Continuous measurement-while-drilling surveying

ABSTRACT

A continuous measurement-while-drilling surveying apparatus includes a fiber optic gyroscope sensitive to rotation about the tool spin axis of the bottom hole assembly and a second fiber optic gyroscope sensitive to rotation of the bottom hole assembly about an axis normal to the tool spin axis. The first gyroscope may be shaped as a torus to accommodate flow of drilling mud through the bottom hole assembly. The outputs of the first and second gyroscopes are processed together with acceleration signals from three accelerometers in a microprocessor which determines the orientation, velocity and position of the bottom hole assembly on a continuous basis.

[0001] This Application claims priority to U.S. Provisional PatentApplication No. 60/262,400 filed Jan. 19, 2001 for “ContinuousMeasurement-While-Drilling Surveying”.

FIELD OF THE INVENTION

[0002] The present invention relates to continuousmeasurement-while-drilling surveying methods and apparatuses.

BACKGROUND OF THE INVENTION

[0003] To obtain hydrocarbons such as oil and gas, boreholes are drilledby rotating a drill bit attached at a drill string end. A largeproportion of the current drilling activity involves directionaldrilling, i.e., drilling deviated and horizontal boreholes to increasethe hydrocarbon production and/or to withdraw additional hydrocarbonsfrom the earth's formations. Modern directional drilling systemsgenerally employ a drill pipe having a drill bit at the bottom that isrotated by a drill motor or a “mud motor”. Pressurized drilling fluid or“mud” or “drilling mud” is pumped into the drill pipe to rotate thedrill motor and to provide lubrication to various members of the drillstring including the drill bit. The drill bit and drill motor form partof what is known as the bottom hole assembly (“BHA”). As required thedrill pipe is rotated by a prime mover, such as a motor, to facilitatedirectional drilling and to drill vertical boreholes.

[0004] Measurement-While-Drilling (MWD) surveying for directional andhorizontal drilling processes is performed to provide the orientationand the position of the BHA [Conti, 1999]. Azimuth, the inclination andthe tool face angles determine the orientation of the BHA, whilelatitude, longitude and altitude determine the position of the BHA. Thealtitude directly determines the true vertical depth of the BHA. Stateof the art MWD surveying techniques are based on magnetic surveyingwhich incorporates three-axis magnetometers and three-axisaccelerometers arranged in three-mutually orthogonal directions. Thethree-axis accelerometers monitor the Earth gravity field to provide theinclination and the tool face angles. This information is combined withthe magnetometer measurements of the Earth magnetic field to provide theazimuth [FIG. 1, Russel et al., 1979].

[0005] The magnetic surveying system determines the BHA orientation atcertain survey stations with the assumption that the error whichmodifies the Earth's magnetic field vector at the surveying instrumentis in the direction of the borehole [Russel et al., 1979]. Thisassumption is justified by installing these magnetometers inside anon-magnetic housing. Such housing system necessitates the use ofnon-magnetic drill collars around the surveying equipment at a costapproaching $30,000 per single installation [Rehm et al., 1989].

[0006] Although simple, the magnetic surveying system suffers fromseveral inadequacies especially within the drilling environment[Thorogood, 1990; Thorogood et al., 1986]. The presence of downhole oredeposits deviate the measurements of the Earth magnetic field even withthe non-magnetic drill collars surrounding the surveying instruments. Inaddition, magnetic surveying tools located in non-magnetic drill collarsare subject to the influence of the other steel components of the drillstring. It has been shown that drill string induced surveying errorincreases with inclination and as the borehole direction approaches theeast-west azimuth. Drill string magnetic interference is particularlynoticeable when inclination exceeds 30°. Although it has been reportedthat the effect of the drill string magnetic interference could bereduced (mitigated but never entirely eliminated) by running longlengths of non-magnetic materials above and below the survey instruments[Grindord et al., 1983], this solution could affect the cost benefits ofthe horizontal drilling technology. The effect of drill string magneticinterference as well as the presence of downhole ore deposits canneither be quantified nor compensated. In addition to these two effects,geomagnetic influences play an important role in the accuracy of themagnetic surveying system [Parkinson, 1983]. Geomagnetic influences aredefined by the variation of the dip angle, the declination and the totalmagnetic field strength with respect to time. The dip angle is the anglebetween the direction of the Earth's magnetic field and the horizontalplane. The declination is the angle between the magnetic north and thetrue north. It was recorded that during any given day at a randomlocation the standard deviations of the dip angle, the declination andthe magnetic field strength were 0.3{square root}, 0.9° and 0.3 μ-Teslarespectively [Parkinson, 1983]. The variation in the geomagnetic fieldis quite significant in relation to the performance capabilities of themagnetic surveying tools currently used. Therefore, geomagnetic effectsmust be taken into account when considering absolute survey accuracy andimpose definite limitations on the accuracy levels that can ultimatelybe achieved.

[0007] Several investigations have been reported to improve the magneticsurveying accuracy. Shiells and Kerridge (2000) introduced theinterpolated in-field referencing method in which absolute localgeomagnetic field data is determined from spot measurements of theEarth's magnetic field. These measurements are taken at local site whichis sufficiently close to the i; drilling site so that the measurementdata is indicative of the Earth's magnetic field at the drilling sitebut is sufficiently remote from the drilling site so that themeasurement data is unaffected by magnetic interference from thedrilling site. This technique compensates for the geomagnetic influencesbut the measurements taken by the magnetometers still suffer from thedrill string magnetic interference. In addition, measurements of thegeomagnetic field in a local site might not be applicable in an areathat has many drilling rigs. Moreover, such technique adds morecomplexity that may give rise to additional cost.

[0008] Hartmann (1998) introduced a method to improve the accuracy ofborehole surveying. This method is based on determining the uncertaintyin the magnetic field measurements by comparing the measured magneticfield with a corresponding known value at specific location. However,this method does not take into account two phenomena. The firstphenomenon is the increase of the drill string magnetic interference athigh inclinations. The second is the presence of ore deposits that canbe found randomly downhole. These two phenomena add additional errorsthat are not taken into consideration when referencing the theoreticalvalue of the magnetic field to the measured value at the known location.

[0009] Dipersio (1995) introduced a new method for the compensation ofgeomagnetic influences. This method depends on matching both thecalculated magnetic field strength and the calculated dip angle to theircorresponding nominal values at a particular geomagnetic location. Thisresults in determining an error-free value of the axial component of themagnetic field which is directly used to determine the azimuth. However,at high inclinations, the drill string magnetic interference generatesadditional errors in the cross-axial directions.

[0010] Several other techniques have been described for magnetic MWDsurveying [Nicholson, 1995; Engebretson, 1992; Helm, 1991; van Dongen etal., 1987; Trowsdale, 1982]. However, it is clear that the majorweaknesses of the present directional sensing instruments stem from theuse of magnetometers for monitoring the azimuth and from the hostileenvironment in which these devices work. The problem encountered withthe use of magnetometers is the presence of massive amount of steelaround the drilling rig. The abundance of ferromagnetic materialnecessitates the separation of the magnetometers by non-magnetic drillcollars. Aside from the cost of utilizing non-magnetic drill collars,their use introduces a second problem. Since the non-magnetic drillcollars impose an additional weight on the drill bit, the surveyingtools are separated from the bearing assembly and the drill bit by about30 feet [Conti, 1989; Rehm at al., 1989]. Elimination of thenon-magnetic drill collars could reduce the distance between theinstrument package and the drill bit by approximately 50%. The thirdproblem associated with the use of magnetometers is their lack ofreliability when used underground due to the deviation of the Earth'smagnetic field from ore deposits.

[0011] It is known to use strapdown inertial navigation systems forborehole surveying using probes. Borehole surveying is performed afterfinishing the drilling process, or sections of the drilling process, forthe purpose of quality assurance. U.S. Pat. No. 4,542,647 describes aborehole inertial guidance system which utilizes two ring lasergyroscopes having their sensitive axes normal to the borehole axis. Inthis case, azimuth is determined by a synthetic rotation signal providedby three-axis accelerometers. While this configuration may be suitablefor borehole mapping, it fails to be useful for MWD surveying. There isconsiderable rotation about the tool-spin axis in a MWD process makingthe use of a synthetic signal to monitor this rotation and thereforeazimuth useless.

[0012] In addition, ring laser gyros are unsuitable for MWD processesbecause they cannot withstand the harsh conditions encountered near theBHA during the drilling process. In particular, the vibration and shockforces associated with the progress of the drill bit impairs theeffectiveness of ring laser gyros.

[0013] Therefore, there is a need in the art for a continuousmeasurement-while-drilling surveying method and apparatus whichmitigates the difficulties of the prior art.

SUMMARY OF THE INVENTION

[0014] The present invention relates to the replacement of the MWDmagnetic surveying system with the technology of inertial navigation.Inertial navigation systems (INS) determine the position, velocity andorientation of a moving body in three-dimensional space by integratingthe measured components of the acceleration and the angular velocity ofthe body [Titteron et al., 1997]. The measurements of acceleration andangular velocity are provided by three-axis accelerometers andthree-axis gyroscopes arranged in three mutually orthogonal directions.

[0015] Currently available 3-gyroscopes/3-accelerometers INS technologycannot provide high accuracy and small dimensions within the same setupand therefore at present, the direct installation of such system insidethe BHA to replace the set of magnetometers is not feasible. Thus, thetraditional high accuracy INS setup is not applicable downhole due tothe limited space inside the BHA and the relatively large dimensions ofthe gyroscopes.

[0016] Therefore, herein disclosed is a method and apparatus fordownhole strapdown INS utilizing only two novel gyroscopes. Onegyroscope has its sensitive axis along the tool spin axis while theother has its sensitive axis normal to the tool spin axis and along theforward direction. This results in the removal of the costlynon-magnetic drill collars because the gyroscopes are insensitive tomagnetic fields. In one embodiment, the two gyroscopes are mountedinside the bearing assembly and are fiber optic gyroscopes (FOG). FOGsare suitable for the downhole drilling application because they do notcontain moving parts, thus providing high reliability with no need forfrequent calibration or maintenance [Kim, 1994; Merhav, 1993; Lefevre,1993]. In addition, the FOGs exhibit low environmental sensitivity sincethey can withstand relatively high temperatures, shocks and vibrations[Noureldin et al., 1999; Noureldin et al., 2000]. Moreover, currentlyavailable FOGs are of small size, with a drift rate less than 1°/hr,long mean time between failure (60,000 hrs), no gravitational effectsand excellent immunity to vibration and shock forces.

[0017] In one embodiment, the FOG that has its sensitive axis along thetool spin axis is prepared as torus to accommodate the flow of drillingfluid through the bearing assembly. The 3 accelerometers—2 FOG system iscapable of providing continuous MWD surveying during the whole drillingprocess, up until the horizontal section of the well, instead of MWDstation-based surveying. In addition, continuous MWD surveying ispossible through an initial radical section of the well with 3accelerometers and a single FOG which has its sensitive axis alignedwith the tool spin axis. The initial radical section of the well may bethat portion with an inclination less than about 45°, preferably lessthan about 30° and more preferably less than about 20°.

[0018] Therefore, in one aspect, the invention comprises a continuousmeasurement-while-drilling surveying apparatus for surveying thedrilling progress of a bottom hole assembly (“BHA”) comprising:

[0019] (a) a housing adapted to be part of the BHA, said housing havinga tool-spin axis;

[0020] (b) a first fiber-optic gyroscope mounted within the housing forgenerating a first angular rotation signal representing angular rotationof the BHA about the tool-spin axis;

[0021] (c) a second fiber-optic gyroscope mounted within the housing forgenerating a second angular rotation signal representing angularrotation of the BHA about an axis normal to the tool-spin axis;

[0022] (d) accelerometer means for generating three acceleration signalsrepresenting the components of acceleration of the BHA along threemutually orthogonal axes;

[0023] (e) first processing means responsive to the acceleration signalsfor determining the angle of the BHA away from the vertical and forgenerating a third angular rotation signal representing rotation of theBHA about an axis normal to the sensitive axes of the first and secondgyroscopes;

[0024] (f) second processing means responsive to the first, second andthird angular rotation signals and the acceleration signals fortransforming signals representing movement of the BHA in a BHAcoordinate system to a earth local-level coordinate system;

[0025] (g) third processing means operatively connected to the secondprocessing means for determining the orientation of the BHA, determiningthe velocity changes of the BHA, updating the velocity components of theBHA and updating the position components of the BHA.

[0026] The first gyroscope may e toroidal. The first, second and thirdprocessing means may be a general purpose computer programmed withappropriate software, firmware, a microcontroller, a microprocessor or aplurality of microprocessors, a digital signal processor or otherhardware or combination of hardware and software known to those skilledin the art.

[0027] In another aspect of the invention, the invention is a method ofcontinuous MWD surveying of a wellbore which includes, on completion, avertical section, a radical section and a horizontal section, using aBHA comprising a first gyroscope having its sensitive axis aligned withthe tool spin axis, a second gyroscope having its sensitive axis on anaxis normal to the tool spin axis and three accelerometers arranged inthree mutually orthogonal directions, said method comprising the stepsof:

[0028] (a) receiving the rotation signals from the first gyroscope andthe second gyroscope and the acceleration signals from the threeaccelerometers;

[0029] (b) establishing the desired azimuth direction;

[0030] (c) determining the pitch angle from the accelerometer signals;

[0031] (d) determining the rotation rate about an axis normal to thesensitive axes of the two gyroscopes by determining the time rate ofchange of the pitch angle;

[0032] (e) determining the effect of the earth rotation and local-level(l-frame) change of orientation to the angular changes monitored alongthe three axes and removing these two effects from the rotation ratesignals;

[0033] (f) transforming the rotation signals from a BHA coordinate frameto the earth local-level coordinate frame;

[0034] (g) determining the pitch, roll and azimuth of the BHA;

[0035] (h) determining the time rate of change of the velocity of theBHA expressed in the l-frame and update the velocity components of theBHA; and

[0036] (i) updating the position components of the BHA.

[0037] In one embodiment, this method is used for an initial portion ofthe radical section of the wellbore. The initial portion of the radicalsection of the well may be that portion with an inclination less thanabout 45°, preferably less than about 30° and more preferably less thanabout 20°. After the initial portion, the method is altered by switchingthe axes orientation of the two gyroscopes. Once the wellbore reachesthe horizontal section and rotary mode drilling is required, thesurveying method may then switch to station-based surveying.

BRIEF DESCRIPTION OF THE DRAWINGS

[0038]FIG. 1 is a block diagram illustrating the functionality of aprior art magnetic surveying system.

[0039]FIG. 2 is a sketch representing the inertial frame.

[0040]FIG. 3 is a sketch representing the Earth-fixed frame.

[0041]FIG. 4 is a sketch representing the local level frame.

[0042]FIG. 5 is a sketch showing the body frame axes arranged inside thebearing assembly of a BHA.

[0043]FIG. 6 shows the change of l-frame orientation along the Earth'ssurface.

[0044]FIG. 7 shows the arrangement of the fiber optic gyroscopes insidea bearing assembly.

[0045]FIG. 8 is a bottom view for the bearing assembly showing the flowof mud through the toroidal FOG.

[0046]FIG. 9 shows the initial arrangement of the FOGs and theaccelerometers with the body frame axes inside the bearing assembly.

[0047]FIG. 10 shows the setup while surveying the radical section until20° inclination.

[0048]FIG. 11 shows the setup after achieving 20° inclination.

[0049]FIG. 12 shows the new axis orientation after 20° inclination.

[0050]FIG. 13 shows the relationship between the pitch and theinclination angles before and after changing the axes orientation.

[0051]FIG. 14 is a block diagram describing the functionality of thesurveying procedure until 20° inclination.

[0052]FIG. 14A is a block diagram describing the functionality of thesurveying procedure until 20° inclination using only the toroidal FOG.

[0053]FIG. 15 is a block diagram describing the functionality of thesurveying procedure after 20° inclination.

[0054]FIG. 16 shows the arrangement of the surveying instruments insidethe bottom hole assembly while drilling the horizontal section of thewell.

[0055]FIG. 17 is a block diagram showing the station-based surveyingprocedure.

[0056]FIG. 18 is a block diagram showing the technique of aided inertialnavigation.

DESCRIPTION OF THE INVENTION

[0057] The present invention provides for a continuousmeasurement-while-drilling surveying method and apparatus. All terms notdefined herein have their common art-recognized meanings.

[0058] Strapdown Intertial Navigation

[0059] There are two common INS categories [Titteron et al., 1997;Britting, 1971]. The first category is the gimbaled INS that maintains astationary platform for the accelerometers by torque motors. Thesesystems use the angular rate measurements of the gyroscopes attached tothe platform to control the motors. The second category is the strapdownINS (SINS) that mathematically transforms the output of theaccelerometers attached to the body into a locally-leveled coordinatesystem before performing integration. These systems use the output ofthe gyroscopes attached to the body to continuously convert from bodycoordinate to locally leveled coordinate.

[0060] Three reference frames are commonly used in the technology ofinertial navigation [Schwarz, 1999; Salychev, 1998]. These are theinertial frame, the earth-fixed frame and the local level frame. Theinertial frame (i-frame) is a non-rotating and non-accelerating frame.The origin of the inertial frame is at the center of mass of the Earth,its Z-axis is parallel to the spin axis of the Earth, its X-axis pointstoward the mean equinoctial cloure and its Y-axis completes a righthanded orthogonal frame (see FIG. 2). The Earth-fixed frame (e-frame)has its origin at the center of mass of the Earth with its Z-axisparallel to the spin axis of the Earth, X-axis pointing towards the meanmeridian of Greenwich and Y-axis completing a right handed orthogonalframe (see FIG. 3). The local level frame (l-frame) has its origin atthe origin of the sensor frame with its Z-axis orthogonal to thereference ellipsoid of the Earth and pointing upward, Y-axis pointingtoward the geodetic north and X-axis completing a right handedorthogonal frame (see FIG. 4). When the inertial sensors (gyroscopes andaccelerometers) are installed inside a moving object, they are orientedin three mutually orthogonal directions known as the body frame(b-frame). The axes of the b-frame point toward the forward (Y),transverse (X) and vertical (Z) directions (see FIG. 5).

[0061] The measurements of linear acceleration and angular velocitiesare taken in the b-frame. These measurements should be transformed intothe corresponding reference frame within which they are processed toprovide the position vector, the velocity vector and the attitudecomponents. The transformation from the b-frame to the l-frame isperformed with the rotation matrix R_(b) ^(l), which is expressed as afunction of the attitude components (azimuth ψ, roll (tool face) φ andpitch (inclination) θ ) as follows: $\begin{matrix}{R_{b}^{l} = \begin{bmatrix}{{\cos \quad {\psi cos\varphi}} - {\sin \quad {\psi sin}\quad {\theta sin\varphi}}} & {{- \sin}\quad {\psi cos\theta}} & {{\cos \quad {\psi sin\theta}} + {\sin \quad {\psi sin\theta cos\varphi}}} \\{{\sin \quad {\psi cos\varphi}} + {\cos \quad {\psi sin\theta sin\varphi}}} & {\cos \quad {\psi cos\theta}} & {{\sin \quad {\psi sin\varphi}} - {\cos \quad {\psi sin\theta cos\varphi}}} \\{{- \cos}\quad {\theta sin\varphi}} & {\sin \quad \theta} & {\cos \quad {\theta cos\varphi}}\end{bmatrix}} & (1)\end{matrix}$

[0062] The position of the BHA as described in the l-frame is expressedin terms of the curvilinear coordinates (latitude φ, longitude λ andaltitude h). The velocity is alternatively expressed by three componentsalong the East direction (V^(e)), north direction (V^(n)) and verticaldirection (V^(u)). Therefore, the time rate of change of the positioncomponents is related to the velocity components as follows [Schwarz etal., 1990]: $\begin{matrix}{{\overset{.}{r}}^{l} = {\begin{pmatrix}\overset{.}{\phi} \\\overset{.}{\lambda} \\\overset{.}{h}\end{pmatrix} = {{\begin{pmatrix}0 & \frac{1}{M + h} & 0 \\\frac{1}{\left( {N + h} \right)\cos \quad \phi} & 0 & 0 \\0 & 0 & 1\end{pmatrix}\begin{pmatrix}V^{e} \\V^{n} \\V^{u}\end{pmatrix}} = {D^{- 1}V^{l}}}}} & (2)\end{matrix}$

[0063] where

[0064] N:is the prime vertical radius of curvature of the Earth'ssurface (East-West);

[0065] M:is the corresponding meridian radius of curvature(North-South);

[0066]_(D) ⁻¹:is a 3-by-3 matrix defined as $\begin{pmatrix}0 & \frac{1}{M + h} & 0 \\\frac{1}{\left( {N + h} \right)\cos \quad \phi} & 0 & 0 \\0 & 0 & 1\end{pmatrix}.$

[0067] The acceleration of the bearing assembly is measured in threemutually orthogonal directions in the b-framef^(b)=(f_(x)f_(y)f_(z))^(T). These measurements are transformed into thel-frame using the rotation matrix R_(b) ^(l) (f^(l)=R_(b) ^(l)f^(b)).However, the acceleration components expressed in the l-frame, f^(l),can not directly give the velocity components of the bearing assemblydue to three reasons. The first reason is due to the rotation rate ofthe Earth about its spin axis (ω^(e)=15°/hr) which is interpreted in thel-frame as ω_(ie) ^(l)=(0ω^(e)cos φ ω^(e)sin φ)^(T). The second reasonis the change of orientation of the l-frame with respect to the Earth asthe bearing assembly penetrates the downhole formation. The l-framechange of orientation along one of the meridians can be shown on FIG. 6.This change of orientation is due to the definition of the local Northand vertical directions. The North direction (N) is tangent to themeridian all the time while the vertical direction (V^(le)) is normal tothe Earth's surface. This effect is interpreted as rotation rate vectorexpressed as ω_(el) ^(l)=(−{dot over (φ)} {dot over (λ)} cos φ{dot over(λ)} sin φ)^(T) which can be rewritten as$\omega_{el}^{l} = \left( {\frac{- V^{n}}{M + h}\frac{V^{e}}{N + h}\frac{V^{e}\tan \quad \phi}{N + h}} \right)^{T}$

[0068] (see equation 2). The third reason is the gravity field of theEarth which is given as g^(l=()0 0−γ)^(T), where γ is obtained from thewell-known normal gravity model given as follows [Schwarz et al., 1999]:

γ=a ₁(1+a ₂sin² φ+a ₃sin⁴φ)+(a ₄ +a ₅sin²φ)h+a ₆ h ²  (3)

[0069] where

[0070] a₁=9.7803267715 m/sec²

[0071] a₂=0.0052790414

[0072] a₃=0.0000232718

[0073] a₄=−0.000003087691089 1/sec²

[0074] a₅=0.000000004397731 1/sec²

[0075] a₆=0.000000000000721 1/(m sec²)

[0076] Then, the time rate of change of the velocity components V^(l) ofthe bearing assembly while penetrating downhole can be expressed asfollows [Schwarz et al., 1999]:

{dot over (V)}^(l) =R _(b) ^(l) f ^(b)−(2Ω_(ie) ^(l)+Ω_(el) ^(l))V ^(l)+g ^(l)  (4)

[0077] where Ω_(ie) ^(l) and Ω_(el) ^(l) are the skew-symmetric matricescorresponding to ω_(ie) ^(l) and ω_(el) ^(l) respectively and they areexpressed as follows: $\begin{matrix}{\Omega_{ie}^{l} = \begin{pmatrix}0 & {{- \omega^{e}}\sin \quad \phi} & {\omega^{e}\cos \quad \phi} \\{\omega^{e}\sin \quad \phi} & 0 & 0 \\{{- \omega^{e}}\cos \quad \phi} & 0 & 0\end{pmatrix}} & (5) \\{\Omega_{el}^{l} = \begin{pmatrix}0 & \frac{{- V^{e}}\tan \quad \phi}{N + h} & \frac{V^{e}}{N + h} \\\frac{V^{e}\tan \quad \phi}{N + h} & 0 & \frac{V^{n}}{M + h} \\\frac{- V^{e}}{N + h} & \frac{- V^{n}}{M + h} & 0\end{pmatrix}} & (6)\end{matrix}$

[0078] The rotation matrix R_(b) ^(l) can be obtained by solving thefollowing differential equations [Schwarz et al., 1999]:

{dot over (R)} _(b) ^(l) =R _(b) ^(l)Ω_(lb) ^(b) =R _(b) ^(l)(Ω_(ib)^(b)−Ω_(il) ^(b))  (7)

[0079] where Ω_(ib) ^(b) is the skew-symmetric matrix of themeasurements of angular velocities provided by the gyroscopes. Ω_(ib)^(b) corresponds to the angular velocity vector ω_(ib)^(b)=(ω_(x)ω_(y)ω_(z))^(T) and is given as: $\begin{matrix}{\Omega_{ib}^{b} = \begin{pmatrix}0 & {- \omega_{z}} & \omega_{y} \\\omega_{z} & 0 & {- \omega_{x}} \\{- \omega_{y}} & \omega_{x} & 0\end{pmatrix}} & (8)\end{matrix}$

[0080] Since the gyroscopes measure both the Earth rotation and thechange in orientation of the l-frame in addition to the angularvelocities of the bearing assembly, the term Ω_(il) ^(b) in equation 7is subtracted from Ω_(ib) ^(b) to remove these two effects. The termΩ_(il) ^(b) consists of two parts. The first part is Ω_(ie) ^(b) whichaccounts for the Earth rotation rate and the second part is Ω_(el) ^(b)which accounts for the orientation change of the l-frame. Therefore,Ω_(il) ^(b) can be written as follows:

Ω_(il) ^(b)=Ω_(ie) ^(b)+Ω_(el) ^(b)  (9)

[0081] Since ω_(ie) ^(b)=R_(l) ^(b)ω_(ie) ^(l) and ω_(el) ^(b)=R_(l)^(b)ω_(el) ^(l), then ω_(il) ^(b)=R_(l) ^(b)(ω_(ie) ^(l)+ω_(el) ^(l))and can be given as follows: $\begin{matrix}{\omega_{il}^{b} = {\left\lbrack {\begin{pmatrix}0 \\{\omega^{e}\cos \quad \phi} \\{\omega^{e}\sin \quad \phi}\end{pmatrix} + \begin{pmatrix}\frac{- V^{n}}{M + h} \\\frac{V^{e}}{N + h} \\\frac{V^{e}\tan \quad \phi}{N + h}\end{pmatrix}} \right\rbrack = {R_{l}^{b}\begin{pmatrix}\frac{- V^{n}}{M + h} \\{\frac{V^{e}}{N + h} + {\omega^{e}\cos \quad \phi}} \\{\frac{V^{e}\tan \quad \phi}{N + h} + {\omega^{e}\sin \quad \phi}}\end{pmatrix}}}} & (10)\end{matrix}$

[0082] Consequently, the corresponding skew-symmetric matrix Ω_(il) ^(b)can be determined.

[0083] When equations 2, 4 and 9 are combined together, they form whatis known as the mechanization equations in the l-frame and they areusually given together as follows [Schwarz et al., 1990; Schwarz et al.,1999]: $\begin{matrix}{\begin{pmatrix}{\overset{.}{r}}^{l} \\{\overset{.}{V}}^{l} \\{\overset{.}{R}}_{b}^{l}\end{pmatrix} = \begin{pmatrix}{D^{- 1}V^{l}} \\{{R_{b}^{l}f^{b}} - {\left( {{2\Omega_{ie}^{l}} + \Omega_{el}^{l}} \right)V^{l}} + g^{l}} \\{R_{b}^{l}\left( {\Omega_{ib}^{b} - \Omega_{il}^{b}} \right)}\end{pmatrix}} & (11)\end{matrix}$

[0084] The inputs to these mechanization equations are the gyroscope andaccelerometer measurements while the outputs are the curvilinearcoordinates, three velocity components and three attitude components.The l-frame is chosen to describe the drill bit motion for threereasons. The first is that the axes of the l-frame are aligned towardsthe local East, North and vertical (normal to the Earth's surface andpointing upward). Using this orientation, the pitch, the roll and theazimuth angles describing the drill bit attitude (inclination, tool faceand azimuth) can be obtained directly as outputs of the mechanizationequations solved in the l-frame. The second reason for choosing thel-frame is that the geodetic coordinates expressing the drill bitlatitude, longitude and depth can be obtained directly by solving themechanization equations described in this frame. Finally, thecomputational errors in the navigation parameters on the horizontalplane (North-East plane) are bounded due to Schuller's effect. Thiseffect stipulates that the inertial system errors of the horizontalplane components are coupled together producing what is called aSchuller loop [Mohammed, 1999; Titterton et al., 1997]. Therefore theseerrors oscillate with frequency called the Schuller frequency (1/5000Hz).

[0085] Apparatus

[0086] At present, the FOG size determines the overall size of theinertial measurement unit. As mentioned earlier, present technologicallevel of FOG-based inertial measurement units do not permit directinstallation of 3-FOGs inside the BHA. Therefore, the present inventioncomprises a novel structure of inertial measurement unit with twospecially designed FOGs in order to accommodate space and precisionrequirements. These two FOGs monitor two rotations rate signals alongthe Y and Z directions of the b-frame axes. The third rotation ratesignal may be generated using the time rate of change of the inclinationangle calculated with the three-axis accelerometer measurements. The twoFOGs and the three-axis accelerometers may be mounted anywhere insidethe drill pipe and the technique provided in the following sections canprocess their measurements to determine the surveying data.

[0087] In one embodiment, the MWD inertial surveying sensors are fittedinside the bearing assembly close behind the drill bit as is shown inFIG. 7. The system comprises a first toroidal FOG (10), a second FOG(12), three accelerometers (14, 16, 18), a microprocessor (20), wirelesscommunication transmitter (22) and a power source such as a battery. Themicroprocessor (20) receives the output signals from the two FOGs (10,12) and the three accelerometers (14, 16, 18) and determines theinertial navigation components as described herein. These components arethen transmitted to the surface for review by a drilling engineer.Alternatively, some or all of the processing may take place at thesurface if the output signals of the inertial sensors are transmitted tothe surface. Standard transmission or telemetry apparatus and methodsmay be used, as is well known in the art of MWD.

[0088] The toroidal FOG (10) is installed with its sensitive axis alongthe tool spin axis. It's shape allows the flow of mud through the drillpipe as is clearly shown in FIG. 8. The dimensions of the toroidal FOG(10) are dictated by the dimensions of the bearing assembly. The secondFOG (12) is installed so that its sensitive axis is normal to the toolspin axis and towards the forward direction (Y-axis) Since its sensitiveaxis is normal to the tool spin axis, this FOG is denoted as the normalFOG (12). Both FOGs should preferably have drift rate less than 1deg./hr and angle random walk of less than 0.1 deg./ hr/ {squareroot}{square root over (Hz)} in 100 Hz bandwidth. The bias drift isdefined as the deviation in the measured rotation rate at constanttemperature. The angle random walk is the broadband random noisecomponent at the FOG output that results either from the shot noise orthe thermal noise in the photodetector.

[0089] The three-axis accelerometers (14, 16, 18) are arranged in threemutually orthogonal directions with one of them having its sensitiveaxis parallel to the tool spin axis.

[0090] Standard wireless communication between the system inside thebearing assembly and the surface MWD processor is provided[Skillingstad, 2000]. A power source such as a battery package providesthe necessary power supply to the two FOGs, the accelerometers, themicroprocessor and the wireless communication system.

[0091] Once the vertical section of the well is established, thehorizontal drilling process involves three main tasks:

[0092] 1. Establishing the desired azimuth direction while the drillpipe is still in the vertical direction.

[0093] 2. Drilling the radical section of the well using steering modeof operation.

[0094] 3. Drilling the horizontal section of the well with using rotarymode of operation.

[0095] The Establishment of the Desired Azimuth Direction

[0096] In order to establish accurately the azimuth direction, precisevalue of the initial azimuth should be known. This initial azimuth valueis determined when the whole system is completely stationary. Therefore,the two FOGs monitor the Earth rotation, while the accelerometersmonitor the Earth gravity. The Earth rotation vector monitored by thegyroscopes, ω_(ie) ^(b), and the Earth rotation vector expressed in thel-frame, ω_(ie) ^(l), are related through the rotation matrix R_(l)^(b)=(R_(b) ^(l))^(T) which transforms the l-frame to the b-frame. Thisrelationship can be written as follows: $\begin{matrix}{\omega_{ie}^{b} = {\left( {\omega_{x}\omega_{y}\omega_{z}} \right)^{T} = {{R_{l}^{b}\omega_{ie}^{l}} = {\begin{pmatrix}{{\cos \quad {\psi cos\varphi}} - {\sin \quad {\psi sin\theta sin\varphi}}} & {{\sin \quad {\psi cos\varphi}} + {\cos \quad {\psi sin\theta sin\varphi}}} & {{- \cos}\quad {\theta sin\varphi}} \\{{- \sin}\quad {\psi cos\theta}} & {\cos \quad {\psi cos\theta}} & {\sin \quad \theta} \\{{\cos \quad {\psi sin\varphi}} + {\sin \quad {\psi sin\theta cos\varphi}}} & {{\sin \quad {\psi sin\varphi}} - {\cos \quad {\psi sin\theta cos\varphi}}} & {\cos \quad {\theta cos\varphi}}\end{pmatrix}\begin{pmatrix}0 \\{\omega^{e}\cos \quad \phi} \\{\omega^{e}\sin \quad \phi}\end{pmatrix}}}}} & (12)\end{matrix}$

[0097] The body frame axes are aligned as shown on FIG. 9 with thenormal FOG having its sensitive axis along the Y-axis and the toroidalFOG having its sensitive axis along the Z-axis (the tool spin axis).Using equation 12, the rotation rate measurement provided by the normalFOG, ω_(y), can be given as function of the Earth rotation rate ω^(e) asfollows:

ω_(y)=(cos ψ cos θ)ω^(e)cos φ+(sin θ)ω^(e)sin φ  (13)

[0098] Consequently, the azimuth can be determined as follows:$\begin{matrix}{{\cos \quad \psi} = \frac{\frac{\omega_{y}}{\omega^{e}\cos \quad \phi} - {\sin \quad {\theta tan\phi}}}{\cos \quad \theta}} & (14)\end{matrix}$

[0099] Since the bottom hole assembly is almost in the verticaldirection, the inclination angle (i.e. the pitch angle, θ) is very closeto zero. Therefore, we can simplify equation 14 by considering sin θ≈0and cos θ>1. The initial azimuth can then be given as: $\begin{matrix}{{\cos \quad \psi} = \frac{\omega_{y}}{\omega^{e}\cos \quad \phi}} & (15)\end{matrix}$

[0100] It should be noted that equation 15 has a singularity when thelatitude angle φ equals to±π/2. This situation corresponds to drillingat the Earth poles that can not exist.

[0101] Therefore the initial azimuth can be determined using equation 15by substituting the angular velocity output of the normal FOG, ω_(y),after removing the offset and the constant drift values. The residualerrors in the measurements provided by this FOG, δω_(y), affect thecalculations of cos ψ and consequently the accuracy of the initialazimuth. The error in cos ψ due to these errorsδω_(y) can be obtained bydifferentiating equation 15: $\begin{matrix}{{\delta \quad \cos \quad \psi} = \frac{\delta \quad \omega_{y}}{\omega^{e}\cos \quad \phi}} & (16)\end{matrix}$

[0102] It is crucial to remove the error δ cos ψ from the cos ψ termcalculated from equation 15 before proceeding to calculate thecorresponding value of the azimuth ψ. This is necessary because of theinversely proportional relationship between the error in the azimuth δψand the azimuth angle itself, ψ, which can be obtained from equation 16as follows: $\begin{matrix}{{\delta\psi} = \frac{\left( \frac{\delta \quad \omega_{y}}{\omega^{e}\cos \quad \phi} \right)}{\sin \quad \psi}} & (17)\end{matrix}$

[0103] Since sin_(Ψ≦)1, the term (1/sin Ψ) in equation 17 tends tomagnify the error δω_(y). When ψ is close to π/2, i.e. approaching theEast direction, the term 1/sin ψ becomes closer to unity and has minoreffect on δΨ. On the other hand, when ψ is close to 0, i.e. approachingthe North direction, the term 1/sin Ψ magnifies the random error δψ. Inorder to avoid this problem, the error term δ cos ψ is determined firstusing equation 16 and removed from the values obtained for cos ψ usingequation 15. Once an accurate estimation of cos ψ is determined, theazimuth can be computed independently from the drilling direction.

[0104] Once the initial azimuth of the bottom hole assembly isdetermined, the drill pipe can be rotated to achieve the desired azimuthdirection. Although the normal FOG is responsible for the determinationof the initial value of the azimuth, the toroidal FOG is responsible forthe establishment of the desired azimuth direction.

[0105] With the BHA oriented along the vertical direction, the drillpipe is rotated along the tool spin axis until the desired azimuthdirection is established. Rotation along the tool spin axis ispreferably the only motion, either angular or translational, performedduring this operation. Under these circumstances, the toroidal FOGmonitors the rotation rate along the tool spin axis in addition to thecomponent of Earth rotation along the vertical direction (ω^(e)sin φ).Therefore, the time rate of change of the azimuth, {dot over (ψ)}, canbe written as follows:

{dot over (ψ)}=ω_(z)−ω^(e)sin φ  (18)

[0106] where ω_(z) is the measurement delivered by the toroidal FOG.

[0107] Equation 18 is solved numerically in real-time to providecontinuous monitoring of the azimuth angle while rotating the drillpipe. This first-order differential equation can be solved numericallyusing the Euler method [Yakowitz et al, 1989] after substituting {dotover (ψ)} by (ψ_(k+1)−ψ_(k))/(Δt), where ψ_(k+1) is the azimuth value attime t_(k+1), ψ_(k) is the azimuth value at time t_(k) andΔt=t_(k+1)−t_(k). Then, the azimuth can be continuously monitored usingthe following equation:

ψ_(k+1)=ψ_(k)+(ω_(z)−ω^(e)sin φ)Δt  (19)

[0108] As the drill pipe continuously rotates about its tool spin axis,the azimuth angle is continuously determined from equation 19. Once thedesired azimuth direction is established, the rotation of the drill pipeis stopped and building the radical section of the well is initiated.

[0109] Surveying the Radical Section of the Well

[0110] During the building of the radical section of the well, thebottom hole assembly starts the drilling process from complete verticaldirection and continues toward a complete horizontal direction. At thebeginning, while the BHA is still in the vertical direction, the bodyframe axes are chosen as shown on FIG. 9. The X, Y, and Z axes pointtoward transverse, forward and vertical directions respectively.Therefore, the Z-axis points toward the tool spin axis. With this choiceof axes orientation, we assume that the motion of the BHA is along theY-axis, and sliding along the Z-axis (FIG. 10). Hence, the toroidal FOGmonitors the changes in the azimuth angle ψ while the normal FOGmonitors the changes in the roll angle φ. However, as the wellboreinclination increases, the Z-axis approaches the horizontal and theY-axis approaches the vertical. The Y-axis is no longer in the forwarddirection and the Z-axis is no longer in the vertical direction. Inaddition, the toroidal FOG is no longer in the horizontal plane tosuccessfully monitor the azimuth. Therefore, at a predeterminedinclination such as 20°, a new orientation of the axes may be set tokeep Y along the forward direction and Z along the vertical direction.This can be achieved by virtually choosing new axes system X¹, Y¹, andZ¹, as is illustrated in FIG. 12. Xl points toward the transversedirection exactly like X while Y¹ points toward the forward direction(i.e. along the tool spin axis) and Z¹ points toward the verticaldirection. This axes orientation is maintained until the end of thedrilling process without any physical change in the mounting of the twoFOGs.

[0111] In a preferred embodiment, there are two surveying methodologiesfor the radical section of the well. The first is applied during aninitial portion of the radical section. As used herein, an “initialportion” is that portion of the well where the inclination is less thanabout 45°, preferablyless than about 30° and more preferably less thanabout 20° inclination. The body frame is defined by the X, Y and Z axesas is shown in FIG. 10. During this initial portion of the well, thetoroidal FOG (with its sensitive axis along the Z-axis) is responsiblefor monitoring the changes in the azimuth angle while the normal FOG(with its sensitive axis along the Y-axis) is responsible for monitoringthe changes in the roll angle. Since the toroidal FOG monitors therotation along the tool spin axis, the tool face angle is directlyobtained from the toroidal FOG output.

[0112] After the initial portion, the tasks of the two FOGs areswitched. In addition, the roll and the tool face angles become the sameand are monitored by the toroidal FOG while the normal FOG monitors thechanges in the azimuth angle.

[0113] It should also be noted that the relationship between the pitchangle and the inclination angle is not the same for the two cases duringand after the initial portion of the well. Before 20° inclination thepitch angle is equal to the inclination angle while after 20°inclination there is π/2 radians phase shift between them. This isbecause the pitch angle is defined as the angle between the forwarddirection (Y-axis) of the moving body and the horizontal plane (FIG.13). Since the definition of the Y-axis is changed at 20° inclination,the definition of the pitch angle is changed accordingly and the π/2radians phase shift is created between the inclination and the pitchangle.

[0114] The pitch angle is determined from the three-axis accelerometers.The determination of the inclination angle is based on the assumptionthat the three-axis accelerometers monitor the gravity components inthree mutually orthogonal directions and their outputs can be integratedtogether to determine the pitch. This assumption is valid as long as thepenetration rate through the downhole formation is relatively small. Theaverage angle-building rate is about 4°/100 ft. for long-radiusapplication, 12°/100 ft. for medium-radius applications and 4⁰/ft. forshort-radius applications. In each of these three categories ofhorizontal drilling, the above assumption can be applied and the pitchcan be determined from the measurements of the accelerometers.

[0115] Surveying the Radical Section Before 20⁰ Inclination

[0116] The surveying technique explained in this sub section providescontinuous surveying based on the measurements delivered by thethree-axis accelerometers, the toroidal FOG and the normal FOG. Thetoroidal FOG determines the rotation rate along the Z-axis, ω_(z), whilethe normal FOG determines the rotation rate along the Y-axis, ω_(y), Therotation rate along the X-axis, ω_(x), can be determined from the timerate of change of the pitch angle. Therefore, the process starts withreceiving the accelerometer measurements that can be combined todetermine the pitch, which, until 20° inclination, is the same as theinclination angle. Based on the assumption of very small penetrationrate through the downhole formation, we assumed that the threeaccelerometers are basically affected by the gravity field of the Earthonly. Therefore, the accelerometer measurements taken at the body frame,g^(b), are related to the gravity field of the Earth expressed in thel-frame, g^(l), through the rotation matrix R_(l) ^(b) as follows:$\begin{matrix}{g^{b} = {\left( \quad \begin{matrix}f_{x} \\f_{y} \\f_{z}\end{matrix}\quad \right) = {{R_{l}^{b}g^{l}} = {\left( \quad \begin{matrix}{{\cos \quad {\psi cos}\quad \varphi} - {\sin \quad \psi \quad \sin \quad \theta \quad \sin \quad \varphi}} & {{\sin \quad {\psi cos}\quad \varphi} + {\cos \quad \psi \quad \sin \quad \theta \quad \sin \quad \varphi}} & {{- \cos}\quad \theta \quad \sin \quad \varphi} \\{{- \sin}\quad \psi \quad \cos \quad \theta} & {\cos \quad \psi \quad \cos \quad \theta} & {\sin \quad \theta} \\{{\cos \quad {\psi sin}\quad \varphi} + {\sin \quad \psi \quad \sin \quad \theta \quad \cos \quad \varphi}} & {{\sin \quad {\psi sin}\quad \varphi} - {\cos \quad \psi \quad \sin \quad \theta \quad \cos \quad \varphi}} & {\cos \quad \theta \quad \cos \quad \varphi}\end{matrix}\quad \right)\begin{pmatrix}0 \\0 \\{- \gamma}\end{pmatrix}}}}} & (20)\end{matrix}$

[0117] where

[0118] f_(x), f_(y) and f_(z) are the measurements provided by the threemutually orthogonal accelerometers.

[0119] γ is the value of gravity obtained from the normal gravity modelgiven in equation 3.

[0120] From equation 20, the accelerometer measurements f_(x), f_(y) andf_(z) can be given in terms of γ as follows:

f _(x)=(cos θ sin φ)γ  (21)

f _(y)×−(sin θ)γ  (22)

f _(z)=−(cos θ cos φ)γ  (23)

[0121] where θ is the pitch (inclination) angle and φ is the roll angle.

[0122] It is clear that the pitch angle can be obtained directly fromequation 22 as long as γ is determined from the normal gravity model(equation 3). The determination of γ depends on the availableinformation for the latitude φ and the altitude h. However, since thenormal gravity model is an approximate model, and in order to preventthe propagation of errors, equations 21 through 23 are manipulated toremove the gravity term and to directly give the pitch angle. Fromequations 21 and 23, we can write

{square root}{square root over (f _(x) ² +f _(z) ²)}=(cos θ)γ  (24)

[0123] Then the pitch angle can be obtained by dividing equation 22 byequation 24. $\begin{matrix}{\theta = {\arctan \left( \sqrt{\frac{f_{y}^{2}}{f_{x}^{2} + f_{z}^{2}}} \right)}} & (25)\end{matrix}$

[0124] Equation 25 gives the pitch angle directly in a way that iscompletely independent of the gravity field of the Earth. Then therotation rate along the X-axis, ω_(x), is determined using the time rateof change of the pitch angle as follows: $\begin{matrix}{\omega_{x} = \frac{\theta_{k + 1} - \theta_{k}}{t_{k + 1} - t_{k}}} & (26)\end{matrix}$

[0125] where θ_(k) and θ_(k+1) are the values of the pitch angle at timet_(k) and t_(k+1) respectively.

[0126] The surveying procedure delivers the navigatio n data (theposition and the orientation) for the bottom hole assembly by solvingthe set of first order differential equations given in equation 11. Thefirst step in solving these equations is the parameterization of therotation matrix R_(b) ^(l) [Schwarz et al., 1999; Salychev, 1998;Titterton et al., 1997]. One method for this purpose is the quatermionapproach. According to Euler's Theorem [Salychev, 1998], the rotation ofa rigid body (represented by the b-frame) with respect to a referenceframe (represented by the l-frame) can be expressed in terms of therotation angle Θ about a fixed axis and the direction cosine of therotation axis that defines the rotation direction. Thus, quaternionparameters Q=(q₁ q₂ q₃ q₄)^(T) are introduced to describe the rotationof the b-frame with respect to the l-frame and they are expressed asfollow: $\begin{matrix}{Q = {\begin{pmatrix}q_{1} \\q_{2} \\q_{3} \\q_{4}\end{pmatrix} = \begin{pmatrix}{\left( {\Theta_{x}/\Theta} \right){\sin \left( {\Theta/2} \right)}} \\{\left( {\Theta_{y}/\Theta} \right){\sin \left( {\Theta/2} \right)}} \\{\left( {\Theta_{z}/\Theta} \right){\sin \left( {\Theta/2} \right)}} \\{\cos \quad \left( {\Theta/2} \right)}\end{pmatrix}}} & (27)\end{matrix}$

[0127] where Θ={square root}{square root over (Θ_(x) ²+Θ_(y) ²+Θ_(z) ²)} is the rotation angle and$\frac{\Theta_{x}}{\Theta},\frac{\Theta_{y}}{\Theta},\frac{\Theta_{z}}{\Theta}$

[0128] are the three direction cosines of the rotation axis with respectto the l-frame.

[0129] The definition of the quaternion parameters described in equation(27) implies that the four quaternion components are not independentbecause q₁ ²+q₂ ²+q₃ ²+q₄ ²=1. This means that only three independentquaternion components are sufficient to describe the rigid bodyrotation. The time rate of change of the quaternion is described by thefollowing first-order differential equation: $\begin{matrix}{\overset{.}{Q} = {\frac{1}{2}{\Omega (\omega)}Q}} & (28)\end{matrix}$

[0130] where Ω(ω) is a skew-symmetric matrix given as follows:$\begin{matrix}{{\Omega (\omega)} = \left( \quad \begin{matrix}0 & \omega_{z} & {- \omega_{y}} & \omega_{x} \\{- \omega_{z}} & 0 & \omega_{x} & \omega_{y} \\\omega_{y} & {- \omega_{x}} & 0 & \omega_{z} \\{- \omega_{x}} & {- \omega_{y}} & {- \omega_{z}} & 0\end{matrix}\quad \right)} & (29)\end{matrix}$

[0131] where ω_(x), ω_(y),ω_(z) are the angular velocities of bodyrotation which are determined by equation 26 for ω_(x) and are monitoredby the normal FOG and the toroidal FOG for ω_(y) and ω_(z) respectively,after compensating for the Earth rotation and the l-frame change oforientation.

[0132] To solve the set of first order differential equations given inequation 28, Euler method can be used to determine the quaternionparameters Q_(k+1) at time t_(k+1) based on the values of the quaternionparameters Q_(k) at time t_(k) as follows: $\begin{matrix}{Q_{k + 1} = {Q_{k} + {\left( {\frac{1}{2}{\Omega \left( \omega_{k} \right)}Q_{k}} \right)\Delta \quad t}}} & (30)\end{matrix}$

[0133] where Δt=t_(k+1)−t_(k).

[0134] Once the quaternion parameters are determined at certain time,the rotation matrix R_(b) ^(l) can be obtained using the followingdirect relationship: $\begin{matrix}\begin{matrix}{R_{b}^{} = \quad {\left( \quad \begin{matrix}R_{11} & R_{12} & R_{13} \\R_{21} & R_{22} & R_{23} \\R_{31} & R_{32} & R_{33}\end{matrix}\quad \right) =}} \\{\quad \left( \quad \begin{matrix}{q_{1}^{2} - q_{2}^{2} - q_{3}^{2} + q_{4}^{2}} & {2\left( {{q_{1}q_{2}} - {q_{3}q_{4}}} \right)} & {2\left( {{q_{1}q_{3}} + {q_{2}q_{4}}} \right)} \\{2\left( {{q_{1}q_{2}} + {q_{3}q_{4}}} \right)} & {{- q_{1}^{2}} + q_{2}^{2} - q_{3}^{2} + q_{4}^{2}} & {2\left( {{q_{2}q_{3}} - {q_{1}q_{4}}} \right)} \\{2\left( {{q_{1}q_{3}} - {q_{2}q_{4}}} \right)} & {2\left( {{q_{2}q_{3}} + {q_{1}q_{4}}} \right)} & {{- q_{1}^{2}} - q_{2}^{2} + q_{3}^{2} + q_{4}^{2}}\end{matrix}\quad \right)}\end{matrix} & (31)\end{matrix}$

[0135] On the other hand, the quaternion parameters can be obtained fromthe rotation matrix R_(b) ^(l) as follows: $\begin{matrix}{\left( \quad \begin{matrix}q_{1} \\q_{2} \\q_{3} \\q_{4}\end{matrix}\quad \right) = \left( \quad \begin{matrix}{0.25{\left( {R_{32} - R_{23}} \right)/q_{4}}} \\{0.25{\left( {R_{13} - R_{31}} \right)/q_{4}}} \\{0.25{\left( {R_{21} - R_{12}} \right)/q_{4}}} \\{0.5\sqrt{1 + R_{11} + R_{22} + R_{33}}}\end{matrix}\quad \right)} & (32)\end{matrix}$

[0136] The quaternion parameters are introduced for the parameterizationof the rotation matrix R_(b) ^(l) for three reasons. The first reason isthat only four differential equations are solved numerically instead ofsix differential equations if we manipulate the rotation matrix R_(b)^(l) directly. The second reason is that the quaternion solution avoidsthe singularity problem that might exist with some other solutionmethods. The third reason is the computational simplicity introduced bythe quaternion.

[0137] One embodiment of the method may be described as follows:

[0138] 1. Receive the measurements from the two FOGs and the threeaccelerometers. Some modern FOGs deliver angular increments at theiroutput instead of rotation rate. Similarly, the accelerometers delivervelocity increments instead of specific force (acceleration)measurements. The angular changes provided by the two FOGs should becompensated for the drift rate d as follows:

Δθ_(y)=Δ{circumflex over (θ_(y))}−d _(y) Δt

Δθ_(z)=Δ{circumflex over (θ_(z))}−d _(z) Δt  (33)

[0139] where

[0140] Δ{circumflex over (θ_(y))} is the raw measurement of the angularincrements provided by the normal FOG;

[0141] Δθ_(y) is the measurement of the angular increments provided bythe normal FOG after compensating for the drift rate d_(y);

[0142] Δ{circumflex over (θ_(z))} is the raw measurement of the angularincrements provided by the toroidal FOG;

[0143] Δθ_(z) is the measurement of the angular increments provided bythe toroidal FOG after compensating for the drift rate d_(z);

[0144] Δt is the difference between the current time sample t_(k+1) andthe previous time sample t_(k).

[0145] Similarly, the velocity increments provided by the accelerometersshould be compensated for the bias terms b as follows:

Δν_(x)=Δ{circumflex over (ν)}_(x) −b _(x) Δt

Δν_(y)=Δ{circumflex over (ν)}_(y) −b _(y) Δt

Δν_(z)=Δ{circumflex over (ν)}_(z) −b _(z) Δt  (34)

[0146] where

[0147] Δ{circumflex over (ν)}_(x) is the raw measurement of the velocityincrements provided by the accelerometer installed along the X-axis;

[0148] Δ{circumflex over (ν)}_(x) Avx is the measurement of the velocityincrements provided by the accelerometer installed along the X-axisafter compensating for the bias b_(x);

[0149] Δ{circumflex over (ν)}_(y) is the raw measurement of the velocityincrements provided by the accelerometer installed along the Y-axis;

[0150] Δν_(y) is the measurement of the velocity increments provided bythe accelerometer installed along the Y-axis after compensating for thebias b_(y);

[0151] Δ{circumflex over (ν)}_(z) is the raw measurement of the velocityincrements provided by the accelerometer installed along the Z-axis;

[0152] Δν_(z) is the measurement of the velocity increments provided bythe accelerometer installed along the Z-axis after compensating for thebias b_(z);

[0153] Δt is the difference between the current time sample t_(k+1) andthe previous time sample t_(k).

[0154] 2. Determine the pitch angle from the accelerometer measurementsas given in equation 25. This equation is slightly modified to suit thevelocity increment output Δν from the accelerometers. This modificationis based on the direct relationship between the specific forcemeasurement f, representing the acceleration, and the velocityincrements Δν given as f=Δν/Δt. Thus, the pitch angle, θ, is determinedas follows: $\begin{matrix}{\theta = {\arctan \left( \sqrt{\frac{\Delta \quad v_{y}^{2}}{{\Delta \quad v_{x}^{2}} + {\Delta \quad v_{z}^{2}}}} \right)}} & (35)\end{matrix}$

[0155] Then, the angular changes along X-direction, Δθ_(x), can beobtained from equation 26 as follows

Δθ_(x)=ω_(x) Δt=[θ(t _(k+1))−θ(t _(k))]  (36)

[0156] 3. Determine the effect of the Earth rotation and the l-framechange of orientation in the angular changes monitored along the threeaxes. The term responsible for the Earth rotation and the l-frame changeof orientation is ω_(il) ^(b) and is given as follows (see equation 10):$\begin{matrix}{\omega_{i}^{b} = {{R_{}^{b}\left( t_{k} \right)}\left( \quad \begin{matrix}\frac{- {V^{n}\left( t_{k} \right)}}{M + h} \\{\frac{V^{e}\left( t_{k} \right)}{N + h} + {\omega^{e}\cos \quad \phi}} \\{\frac{{V^{e}\left( t_{k} \right)}\tan \quad \phi}{N + h} + {\omega^{e}\sin \quad \phi}}\end{matrix}\quad \right)}} & (37)\end{matrix}$

[0157] Then the angular changes corresponding to ω_(il) ^(b) can bedetermined as:

θ_(il) ^(b)=ω_(il) ^(b) Δt  (38)

[0158] Consequently, the vector of measured angular increments θ_(ib)^(b)=(Δθ_(x) Δθ_(y) Δθ_(z))^(T) is compensated for θ_(il) ^(b) todetermine the actual bottom hole assembly angular changes θ_(ib) ^(b) asfollows:

θ_(lb) ^(b)=θ_(ib) ^(b)−θ_(il) ^(b)=(δθ_(x) δθ_(y) δθ_(z))^(T)  (39)

[0159] 4. Update the quaternion. The time rate of change of thequaternion as given in equation 28 corresponds to the time rate ofchange of the rotation matrix R_(b) ^(l) as given in equation 7. Thequaternion update follows equation 30. The time interval Δt ismultiplied by the elements of the skew-symmetric matrix of angularvelocities Ω, so that ω_(x)Δt=δθ_(x), ω_(y)Δt=δθ_(y) and ω_(z)Δt=δθ_(z).The quaternion update equation is then given as follows: $\begin{matrix}\begin{matrix}{\left( \quad \begin{matrix}{q_{1}\left( t_{k + 1} \right)} \\{q_{2}\left( t_{k + 1} \right)} \\{q_{3}\left( t_{k + 1} \right)} \\{q_{4}\left( t_{k + 1} \right)}\end{matrix}\quad \right) = \quad {\left( \quad \begin{matrix}{q_{1}\left( t_{k} \right)} \\{q_{2}\left( t_{k} \right)} \\{q_{3}\left( t_{k} \right)} \\{q_{4}\left( t_{k} \right)}\end{matrix}\quad \right) +}} \\{\quad {\frac{1}{2}\left( \quad \begin{matrix}0 & {\delta \quad \theta_{z}} & {{- \delta}\quad \theta_{y}} & {\delta \quad \theta_{x}} \\{{- \delta}\quad \theta_{z}} & 0 & {\delta \quad \theta_{x}} & {\delta \quad \theta_{y}} \\{\delta \quad \theta_{y}} & {{- \delta}\quad \theta_{x}} & 0 & {\delta \quad \theta_{z}} \\{{- \delta}\quad \theta_{x}} & {{- \delta}\quad \theta_{y}} & {{- \delta}\quad \theta_{z}} & 0\end{matrix}\quad \right)\left( \quad \begin{matrix}{q_{1}\left( t_{k} \right)} \\{q_{2}\left( t_{k} \right)} \\{q_{3}\left( t_{k} \right)} \\{q_{4}\left( t_{k} \right)}\end{matrix}\quad \right)}}\end{matrix} & (40)\end{matrix}$

[0160] where the initial value of the quaternion Q(t₀)=(q₁(t₀) q₂(t₀)q₃(t₀) q₄(t₀) is chosen to be equal to (0 0 0 1)^(T).

[0161] Consequently, the rotation matrix R_(b) ^(l) can be determinedfrom its direct relationship to the quaternion Q =(q₁ q₂ q₃ q₄)^(T) asgiven by equation 31. From the Re definition given in equation 1, thepitch (inclination) is determined as θ=arcs in(R₃₂), the roll isdetermined as φ arctan(−R₃₁/R₃₃) and the azimuth (tool face) isdetermined as ψ=arctan(−R₁₂/R₂₂)

[0162] 5. Determine the changes in velocity of the bottom hole assemblyusing equation 4 that describes the time rate of change of the velocityof the bottom hole assembly expressed in the l-frame. This equation canbe rewritten as follows:$\frac{\Delta \quad V^{}}{\Delta \quad t} = {{R_{b}^{}f^{b}} - {\left( {{2\Omega_{i\quad e}^{}} + \Omega_{e\quad }^{}} \right)V^{}} + g^{}}$

[0163] or

ΔV ^(l) =R _(b) ^(l) f ^(b) Δt−(2Ω_(ie) ^(l)+Ω_(el) ^(l))V ^(l) Δt+g^(l) Δt  (41)

[0164] Since the specific force measurement f^(b)=(f_(x) f_(y)f_(z))^(T) is related to the velocity increment measurementsΔν^(b)=(Δν_(x) Δν_(y) Δν_(z))^(T) with f^(b)=Δν^(b)Δt, equation 41 canbe given as:

ΔV ^(l)(t _(k+1))=R _(b) ^(l)Δν^(b)(t _(k+1))−(2Ω_(ie) ^(l)+Ω_(el)^(l))V ^(l)(t _(k))Δt+g ^(l) Δt  (42)

[0165] where Ω_(ie) ^(l) and Ω_(el) ^(l) are given in equations 5 and 6respectively, while g^(l) is the gravity field vector and is defined asshown in equation 3.

[0166] The first term on the right-hand side of equation 42 is themeasured velocity increments transformed to the l-frame. The second termis the Coriolis correction that compensates for the Earth rotation andthe l-frame change of orientation. The third tem is the gravitycorrection.

[0167] The velocity components at the current time epoch t_(k+1) arecomputed using the modified Euler formula as follows: $\begin{matrix}{{V^{}\left( t_{k + 1} \right)} = {{V^{}\left( t_{k} \right)} + {\frac{1}{2}\left( {{\Delta \quad {V^{}\left( t_{k} \right)}} + {\Delta \quad {V^{}\left( t_{k + 1} \right)}}} \right)}}} & (43)\end{matrix}$

[0168] where V^(l)=(V^(e) V^(n) V^(u)) are the velocity components asdescribed in the l-frame.

[0169] 6. Determine the altitude of the bottom hole assembly, h, (i.e.the true vertical depth TVD) using the vertical component of thevelocity vector, V^(u). Since h=V^(u) according to equation 2, the TVDcan be computed at time instant t_(k+1) using the modified Euler formulaas follows: $\begin{matrix}{{h\left( t_{k + 1} \right)} = {{h\left( t_{k} \right)} + {\frac{1}{2}\left( {{V^{u}\left( t_{k + 1} \right)} + {V^{u}\left( t_{k} \right)}} \right)\Delta \quad t}}} & (44)\end{matrix}$

[0170] Once the TVD is determined, the latitude φ, and the longitude λ,can be computed. According to equation 2, φ and λ can be computed attime instant (t_(k+1)) using the modified Euler formula as follows:$\begin{matrix}{{\phi \left( t_{k + 1} \right)} = {{\phi \left( t_{k} \right)} + {\frac{1}{2}\frac{\left( {{V^{n}\left( t_{k + 1} \right)} + {V^{n}\left( t_{k} \right)}} \right)}{M + h}\Delta \quad t}}} & (45) \\{{\lambda \quad \left( t_{k + 1} \right)} = {{\lambda \quad \left( t_{k} \right)} + {\frac{1}{2}\frac{\left( {{V^{e}\left( t_{k + 1} \right)} + {V^{e}\left( t_{k} \right)}} \right)}{\left( {N + h} \right)\cos \quad \phi}\Delta \quad t}}} & (46)\end{matrix}$

[0171] It should be noted that the computation of the longitude, λ, willsuffer from singularities if this procedure is performed closer to thepoles where the latitude φ becomes equal to π/2.

[0172] The six steps described above are repeated for each new set ofmeasurements received from the two FOGs and the three accelerometersuntil 20° inclination. A block diagram describing the functionality ofthis procedure is shown on FIG. 14.

[0173] In an alternative embodiment of the invention, surveying may beaccomplished with only the toroidal FOG and the three accelerometersuntil 20° inclination. A block diagram describing the functionality ofthis procedure is shown in FIG. 14A. In this case, the toroidal FOG maymonitor both azimuth and roll angles because of the near verticality ofthe BHA. As the pitch angle increases, the unity between azimuth androll decreases and the ability of the toroidal FOG to monitor roll anglediminishes.

[0174] Surveying the Radical Section After 20° Inclination

[0175] As shown earlier on FIG. 13, change of axes orientation may takeplace at 20° inclination in one embodiment of the invention. The Y¹-axisis now oriented along the tool spin axis while Z¹-axis is oriented alongthe orthogonal direction. The X¹-axis remains collinear to the oldX-axis. The initial values for the navigation parameters (positioncomponents, velocity components and attitude components) for surveyingthis part of the radical section of the well can be obtained from thelast values determined just before changing the axes orientation. Itshould be reiterated that the location of the two FOGs and theaccelerometers would not physically change due to the change of axesorientation. Once the procedure for the first part of the radicalsection of the well (before 20° inclination) indicates an inclination of20°, the output of the toroidal FOG is interpreted as Δ{circumflex over(θ)}_(y), since its sensitive axis is reoriented along Y¹-axis.Similarly, the output of the normal FOG is interpreted as Δ{circumflexover (θ)}_(z) since its sensitive axis is reoriented along the Z¹-axis(see FIG. 12). In addition, the output of the accelerometer installedalong the tool spin axis (Y¹-axis) is interpreted as Δ{circumflex over(ν)}_(y) instead of Δ{circumflex over (ν)}_(z) while the outpout of theaccelerometer installed along the orthogonal direction (Z¹-axis) isinterpreted as Δ{circumflex over (ν)}_(z) instead of Δ{circumflex over(ν)}_(y). The third accelerometer installed along the X¹-axis has itsoutput interpreted as Δ{circumflex over (ν)}_(x) the same as in thefirst part of the drilling of the radical section of the well.

[0176] The navigation procedure for this part of the radical section ofthe well is exactly the same as the procedure described in the previoussection except for three changes. The first change is the differentinterpretation of the measurements provided by the two FOGs and thethree accelerometers. The second change is that the tool face anglebecomes no longer equivalent to the azimuth. Instead, it becomesequivalent to the roll. The third change is that the inclination is nolonger equivalent to the pitch (θ) and is given as π/2−θ. The navigationprocedure for surveying after 20° inclination is shown schematically inthe block diagram shown as FIG. 15.

[0177] Surveying the Horizontal Section of the Well

[0178] Surveying the horizontal section of the well is different fromthe radical section due to the rotary mode of drilling during which thewhole drill pipe starts to continuously rotate about its spin axis.Therefore, continuous surveying can no longer be performed andstation-based surveying procedure should be prepared. The station-basedsurveying is based on stopping the drilling operation at certain surveystations. Since the surveying setup is completely stationary at eachstation, the two FOGs are affected by the Earth rotation rate and thethree accelerometers are affected by the Earth gravity field. Thearrangement of the surveying instruments inside the bottom hole assemblyis shown on FIG. 16. The relationship between the Earth rotation ratevector monitored by the FOGs and the same vector defined at the l-frameis given using the rotation matrix R_(l) ^(b) as follows:$\begin{matrix}{\omega_{ie}^{b} = {\begin{pmatrix}\omega_{x} \\\omega_{y} \\\omega_{z}\end{pmatrix} = {{R_{l}^{b}\omega_{ie}^{l}} = {\left( \quad \begin{matrix}{{\cos \quad \psi \quad \cos \quad \varphi} - {\sin \quad \psi \quad \sin \quad \theta \quad \sin \quad \varphi}} & {{\sin \quad \psi \quad \cos \quad \varphi} + {\cos \quad \psi \quad \sin \quad \theta \quad \sin \quad \varphi}} & {{- \cos}\quad \theta \quad \sin \quad \varphi} \\{{- \sin}\quad \psi \quad \cos \quad \theta} & {\cos \quad \psi \quad \cos \quad \theta} & {\sin \quad \theta} \\{{\cos \quad \psi \quad \sin \quad \varphi} + {\sin \quad \psi \quad \sin \quad \theta \quad \cos \quad \varphi}} & {{\sin \quad \psi \quad \sin \quad \varphi} - {\cos \quad \psi \quad \sin \quad \theta \quad \cos \quad \varphi}} & {\cos \quad \theta \quad \cos \quad \varphi}\end{matrix}\quad \right)\begin{pmatrix}0 \\{\omega^{e}\cos \quad \phi} \\{\omega^{e}\sin \quad \phi}\end{pmatrix}}}}} & (47)\end{matrix}$

[0179] The angular velocity ω_(y) is monitored by the toroidal FOG andcan be obtained from the corresponding angular increment output asfollows: $\begin{matrix}{\omega_{y} = \frac{{\Delta\theta}_{y}}{\Delta \quad t}} & (48)\end{matrix}$

[0180] The component of the Earth rotation rate monitored by thetoroidal FOG can be obtained from equation (47) and written as:

ω_(y)=(cos Ψ cos θ)ω^(e) cos φ+(sin θ)ω^(e) sin φ  (49)

[0181] This equation is similar to equation 13 for the establishment ofthe desired azimuth direction. The only difference is that the angularvelocity measurement, ω_(y), is provided by the toroidal FOG instead ofthe normal FOG. With simple mathematical manipulation, equation 49 canbe rewritten as: $\begin{matrix}{{\cos \quad \psi} = \frac{\frac{\omega_{y}}{\omega^{e}\cos \quad \phi} - {\sin \quad \theta \quad \tan \quad \phi}}{\cos \quad \theta}} & (50)\end{matrix}$

[0182] Since the drilling process for this section of the well is alongthe horizontal plane, we can assume that the pitch angle is equal tozero and equation 50 can be simplified. However, general solution ispresented in this section to suit a special case that might exist whiledrilling the radical section of the well. This case involves drillingsome parts of the radical section with rotary mode of drilling insteadof utilizing the steering mode. In this special case the pitch anglecannot be assumed equal to zero.

[0183] To determine the azimuth angle from equation 50, the pitch angle,θ, should be computed. Moreover, complete station based surveyingrequires the determination of the inclination angle (π/2−θ) and the toolface angle (the roll angle, φ). This information is determined bymanipulating the accelerometer measurements. The accelerometermeasurements (f_(x) f_(y) f_(z))^(T) are related to the Earth gravityfiled through the rotation matrix R_(l) ^(b). This relationship can bewritten as: $\begin{matrix}{g^{b} = {\begin{pmatrix}f_{x} \\f_{y} \\f_{z}\end{pmatrix} = {{R_{l}^{b}g^{l}} = {\left( \quad \begin{matrix}{{\cos \quad \psi \quad \cos \quad \varphi} - {\sin \quad \psi \quad \sin \quad \theta \quad \sin \quad \varphi}} & {{\sin \quad \psi \quad \cos \quad \varphi} + {\cos \quad \psi \quad \sin \quad \theta \quad \sin \quad \varphi}} & {{- \cos}\quad \theta \quad \sin \quad \varphi} \\{{- \sin}\quad \psi \quad \cos \quad \theta} & {\cos \quad \psi \quad \cos \quad \theta} & {\sin \quad \theta} \\{{\cos \quad \psi \quad \sin \quad \varphi} + {\sin \quad \psi \quad \sin \quad \theta \quad \cos \quad \varphi}} & {{\sin \quad \psi \quad \sin \quad \varphi} - {\cos \quad \psi \quad \sin \quad \theta \quad \cos \quad \varphi}} & {\cos \quad \theta \quad \cos \quad \varphi}\end{matrix}\quad \right)\begin{pmatrix}0 \\0 \\{- \gamma}\end{pmatrix}}}}} & (51)\end{matrix}$

[0184] The expressions for f_(x), f_(y) and f_(z) can be expressed interms of the pitch and the roll angles by manipulating equation 51 asfollows:

f _(x)=(cos θ sin φ)γ  (52)

f _(y)=−(sin θ)γ  (53)

f _(z)=−(cos θ cos φ)γ  (54)

[0185] From equations (52-54) the pitch and the roll angles can beexpresses as: $\begin{matrix}{\theta = {\arctan \left( \sqrt{\frac{f_{y}^{2}}{f_{x}^{2} + f_{z}^{2}}} \right)}} & (55) \\{\varphi = {\arctan \left( {- \frac{f_{x}}{f_{z}}} \right)}} & (56)\end{matrix}$

[0186] The specific force measurements f_(x), f_(y) and f_(z) areobtained from the velocity increment outputs of the accelerometersΔν_(x), Δv_(y) and Δν_(z) as follows: $\begin{matrix}{{f_{x} = \frac{\Delta \quad v_{x}}{\Delta \quad t}},{f_{y} = \frac{\Delta \quad v_{y}}{\Delta \quad t}},{f_{z} = \frac{\Delta \quad v_{z}}{\Delta \quad t}}} & (57)\end{matrix}$

[0187] Therefore, station-based surveying during the horizontal sectionof the well delivers only the orientation of the bottom hole assemblyusing equations 50, 55 and 56 to provide the azimuth, the pitch and theroll respectively. The inclination is then equal to π/2−θ and the toolface is exactly the roll.

[0188] Error Analysis for the Station-based Surveying

[0189] Although the measurements from the toroidal FOG and theaccelerometers are compensated for all bias terms, there are still someresidual errors that affect the accuracy of the attitude angles. Theseresidual errors include the cross-axis sensitivity,the bias calibrationerror, the bias temperature shift, and the scale factor calibrationerror. Some typical values are assigned to these errors and they aredefined by the manufacturer.

[0190] The error components in the specific force measurements can beobtained by differentiating equations 52, 53 and 54.

δf _(x)=(−sin θ sin φδθ+cos θ cos φδφ)γ  (58)

δf _(y)=−(cos θδθ)γ  (59)

δf _(z)=(sin θ cos φδθ+cos θ sin φδφ)γ  (60)

[0191] The error in the pitch angle, δθ, can be obtained from equation59 and is expressed as follows: $\begin{matrix}{{\delta \quad \theta} = {- \frac{\left( \frac{\delta \quad f_{y}}{\gamma} \right)}{\cos \quad \theta}}} & (61)\end{matrix}$

[0192] Using equations 58 and 60, (δf_(x))²+(δf_(z))² can be given asfollows:

(δf _(x))²+(δf _(z))²=(sin ²θ(δθ)²+cos ²θ(δφ)²)γ²  (62)

[0193] Equation (62) can be manipulated with the substitution of theerror in the pitch angle from equation (61) to give the error in theroll angle δφ as follows: $\begin{matrix}{{\delta \quad \varphi} = \frac{\left( \frac{\sqrt{{\delta \quad f_{x}^{2}} + {\delta \quad f_{z}^{2}} - {\delta \quad f_{y}^{2}\tan^{2}\theta}}}{\gamma} \right)}{\cos \quad \theta}} & (63)\end{matrix}$

[0194] The errors in both the pitch and the roll are determined fromequations 61 and 63 and are used to correct the values of the pitch andthe roll angles obtained by equations 55 and 56. The correct value ofthe pitch angle is then used in equation 50 to determine the azimuth.However, any residual errors existing in the toroidal FOG affect theazimuth accuracy. The azimuth error can be determined from equation 49.If this equation is differentiated with the assumption that the pitchangle is free from error, the toroidal FOG error can be expressed as:

δω_(y)=δ(cos Ψ)cos θ(ω^(e) cos φ)  (64)

[0195] Then the error term δ(cos Ψ) can be written as $\begin{matrix}{{\delta \quad \left( {\cos \quad \psi} \right)} = \frac{\left( \frac{\delta \quad \omega_{y}}{\left( {\omega^{e}\cos \quad \phi} \right)} \right)}{\cos \quad \theta}} & (65)\end{matrix}$

[0196] As discussed above regarding the establishment of the desiredazimuth direction, the error in the cosΨ term should be determined usingequation 65 and completely removed before the value of the azimuth angleis computed. This is to avoid the magnification of the azimuth error ifcalculated directly as follows: $\begin{matrix}{{\delta \quad \psi} = {{- \frac{\left( \frac{\delta \quad \omega_{y}}{\left( {\omega^{e}\cos \quad \phi} \right)} \right)}{\cos \quad \theta}} \cdot \frac{1}{\sin \quad \psi}}} & (66)\end{matrix}$

[0197] It can be inferred from equation 66 that the term 1/sin Ψ largelymagnifies the azimuth error, δΨ, as the bottom hole assembly approachesthe North-South direction (i.e. when Ψ≈0). In addition, singularityproblems exist in the calculation of δΨif δ=0 . Therefore, to avoid allthese problems, the error term δ(cos Ψ) is computed prior to determiningthe azimuth.

[0198] The whole station-based surveying procedure can be shown as ablock diagram in FIG. 17. It should be mentioned that this proceduregives coarse computation of the attitude angles. Precise computationshould be applied to improve the estimation accuracy. The technique forprecise computation is described in the following sections.

[0199] Error Analysis of the Motion of the Bottom Hole Assembly

[0200] In practical systems employing the technology of strapdowninertial navigation the accuracy is limited by the errors originating inthe raw measurements as well as by the imperfections in the variouscomponents which are combined to build the system. These error sourcesare categorized as inertial sensor errors and computational errors[Stephenson et al., 1992; Schwarz et al., 1990]. Any lack of precisionin the measurements of an inertial navigation system is passed from oneestimate to the next with the overall uncertainty in the precision ofthe estimated quantity drifting with time. This explains why theaccuracy of the INS deteriorates in the long term. For MWD surveying,the drilling process continues for several hours. Therefore, errormodels are required for analysis and estimation of the various errorsources associated with the inertial system proposed in the previoussection.

[0201] Errors in dynamic systems are variable in time and thereforedescribed by differential equations. Linearization of the non-lineardynamic system is the most common approach for deriving a set of lineardifferential equations which define the error states of the dynamicsystem. The error states expressed in the t -frame for the mechanizationequations given in the previous section are defined as follows:

e ^(l)=(δφ,δλ,δh,δV _(e) ,δV _(n) ,δV _(u),δθ,δφ,δΨ,δω_(x),δω_(y),δω_(z),δf _(x) ,δf _(y) ,δf _(z))  (67)

[0202] where

[0203] δφ,δλ,δh are the coordinate errors;

[0204] δV^(e),δV^(n),δV^(u) are the velocity errors;

[0205] δθ,δφ,δΨ are the attitude errors;

[0206] δω_(x),δω_(y),δω_(z) are the errors in the angular velocitymeasurements;

[0207] δf_(x),δf_(y),δf_(z) are the errors in the specific forcemeasurements.

[0208] Attitude Errors

[0209] The attitude errors are defined as the misalignment due to theorthogonal transformation between the b-frame and the l-frame. Theseerrors (δθ, δφ and δΨ) depend on two main sources. The first source arethe errors in the measurement of angular velocities provided by theFOGs. The second source are the errors generated due to the Earthrotation and the change of orientation of the l-frame. The expressiondescribing the time rate of change of the attitude errors is given asfollows:

ε^(.l)=−Ω_(il) ^(l)ε^(l)−δω_(il) ^(l) +R _(b) ^(l) d  (68)

[0210] where

[0211] ε^(l)=(δθ δφ δω)^(T) is the attitude error;

[0212] Ω_(il) ^(l)=Ω_(ie) ^(l)+Ω_(el) ^(l) is the skew symmetric matrixfor the rotation of the l-frame with respect to the i-frame;

[0213] δω_(il) ^(l)=δω_(ie) ^(l)+δω_(el) ^(l) is the angular velocityerror vector corresponding to ω_(il) ^(l);

[0214] d=(δω_(x) δω_(y) δω_(z))^(T) is the vector of angular velocitiesmeasurement errors (FOG drift rates).

[0215] The second term in the right hand side of equation 68, δω_(il)^(l), depends implicitly on the velocity errorsδV^(l)=(δV^(e)δV^(n)δV^(u))^(T) and the coordinate errors δr^(l)=(δφ δλδh)^(T). This angular velocity error consists of two terms, δω_(ie) ^(l)and δω_(el) ^(l), which are given as: $\begin{matrix}{{\delta \quad \omega_{ie}^{l}} = \quad {\left( \quad \begin{matrix}0 & 0 & 0 \\{{- \omega^{e}}\sin \quad \phi} & 0 & 0 \\{\omega^{e}\cos \quad \phi} & 0 & 0\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad \phi} \\{\delta \quad \lambda} \\{\delta \quad h}\end{pmatrix}}} & (69) \\{{\delta \quad \omega_{e\quad }^{l}} = {{\left( \quad \begin{matrix}0 & 0 & \frac{V^{n}}{\left( {M + h} \right)^{2}} \\0 & 0 & \frac{- V^{e}}{\left( {N + h} \right)^{2}} \\\frac{V^{e}}{\left( {N + h} \right)\cos^{2}\phi} & 0 & \frac{{- V^{e}}\tan \quad \phi}{\left( {N + h} \right)^{2}}\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad \phi} \\{\delta \quad \lambda} \\{\delta \quad h}\end{pmatrix}} + {\left( \quad \begin{matrix}0 & \frac{- 1}{M + h} & 0 \\\frac{1}{N + h} & 0 & 0 \\\frac{\tan \quad \phi}{N + h} & 0 & 0\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad V^{e}} \\{\delta \quad V^{n}} \\{\delta \quad V^{u}}\end{pmatrix}}}} & (70)\end{matrix}$

[0216] With the substitution of the expressions of δω_(ie) ^(l)(equation 69) and δω_(el) ^(l) (equation 70), the attitude errorequation (equation 68) can be rewritten as: $\begin{matrix}\begin{matrix}{\begin{pmatrix}{\delta \quad \overset{.}{\theta}} \\{\delta \quad \overset{.}{\varphi}} \\{\delta \quad \overset{.}{\psi}}\end{pmatrix} = \quad {{\left( \quad \begin{matrix}0 & {{{\omega \quad}^{e}\sin \quad \phi} + \frac{V^{e}\tan \quad \phi}{N + h}} & {{{- \omega^{e}}\cos \quad \phi} - \frac{V^{e}}{N + h}} \\{{{- \omega^{e}}\sin \quad \phi} - \frac{V^{e}\tan \quad \phi}{N + h}} & 0 & {- \frac{V^{n}}{M + h}} \\{{\omega^{e}\cos \quad \phi} + \frac{V^{e}}{N + h}} & \frac{V^{n}}{M + h} & 0\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad \theta} \\{\delta \quad \varphi} \\{\delta \quad \psi}\end{pmatrix}} +}} \\{\quad {{\left( \quad \begin{matrix}0 & 0 & \frac{- V^{n}}{\left( {M + h} \right)^{2}} \\{\omega^{e}\sin \quad \phi} & 0 & \frac{V^{e}}{\left( {N + h} \right)^{2}} \\{{{- \omega^{e}}\cos \quad \phi} - \frac{V^{e}}{\left( {N + h} \right)\cos^{2}\phi}} & 0 & \frac{V^{e}\tan \quad \phi}{\left( {N + h} \right)^{2}}\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad \phi} \\{\delta \quad \lambda} \\{\delta \quad h}\end{pmatrix}} +}} \\{\quad {{\left( \quad \begin{matrix}0 & \frac{1}{M + h} & 0 \\\frac{- 1}{N + h} & 0 & 0 \\\frac{{- \tan}\quad \phi}{N + h} & 0 & 0\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad V^{e}} \\{\delta \quad V^{n}} \\{\delta \quad V^{u}}\end{pmatrix}} + {R_{b}^{l}d}}}\end{matrix} & (71)\end{matrix}$

[0217] Coordinate Errors

[0218] Since the coordinates of the bottom hole assembly depend only onthe velocity components V^(e), V^(n) and V^(u), the coordinate errorsδφ, δλ and δh are given directly in terms of velocity errors δV^(e),δV^(n) and δV^(u). The expression for the coordinate errors is obtainedby differentiating both sides of equation 2 that describes therelationship between the time rate of change of the coordinatecomponents and the velocity components. Since this relationship is givenby the matrix D⁻¹ with its elements depending on the coordinatecomponents φ and h, the coordinate errors δφ, δλ and δh are written asfollows: $\begin{matrix}\begin{matrix}{{\delta \quad {\overset{.}{r}}^{l}} = \quad {{\left( \quad \begin{matrix}0 & \frac{1}{M + h} & 0 \\\frac{1}{\left( {N + h} \right)\cos \quad \phi} & 0 & 0 \\0 & 0 & 1\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad V^{e}} \\{\delta \quad V^{n}} \\{\delta \quad V^{u}}\end{pmatrix}} +}} \\{\quad {\left( \quad \begin{matrix}0 & 0 & \frac{- V^{n}}{\left( {M + h} \right)^{2}} \\\frac{V^{e}\tan \quad \phi}{\left( {N + h} \right)\cos \quad \phi} & 0 & \frac{- V^{e}}{\left( {N + h} \right)^{2}\cos \quad \phi} \\0 & 0 & 0\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad \phi} \\{\delta \quad \lambda} \\{\delta \quad h}\end{pmatrix}}}\end{matrix} & (72)\end{matrix}$

[0219] Velocity Errors

[0220] The complete implementation of the l-frame mechanizationequations requires the computation of the three velocity componentsV^(e), V^(n) and V^(u) using equation 4. The corresponding errorequations are obtained by differentiating equation 4 and are given as:$\begin{matrix}\begin{matrix}{\begin{pmatrix}{\delta \quad {\overset{.}{V}}^{e}} \\{\delta \quad {\overset{.}{V}}^{n}} \\{\delta \quad {\overset{.}{V}}^{u}}\end{pmatrix} = \quad {{\left( \quad \begin{matrix}{{\left( {2\omega^{e}\sin \quad \phi} \right)V^{u}} + {\left( {2\omega^{e}\cos \quad \phi} \right)V^{n}} + \frac{V^{n}V^{e}}{\left( {N + h} \right)\cos^{2}\phi}} & 0 & 0 \\{{\left( {{- 2}\omega^{e}\cos \quad \phi} \right)V^{e}} - \frac{\left( V^{e} \right)^{2}}{\left( {N + h} \right)\cos^{2}\phi}} & 0 & 0 \\{\left( {{- 2}\omega^{e}\sin \quad \phi} \right)V^{e}} & 0 & \frac{2\quad \gamma}{R}\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad \phi} \\{\delta \quad \lambda} \\{\delta \quad h}\end{pmatrix}} +}} \\{\quad {{\left( \quad \begin{matrix}0 & f_{u} & {- f_{n}} \\{- f_{u}} & 0 & f_{e} \\f_{n} & {- f_{e}} & 0\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad \theta} \\{\delta \quad \varphi} \\{\delta \quad \psi}\end{pmatrix}} +}} \\{\quad {{\left( \quad \begin{matrix}{\frac{- V_{u}}{N + h} + {\frac{V_{n}}{N + h}\tan \quad \phi}} & {{2\omega^{e}\sin \quad \phi} + \frac{V^{e}\tan \quad \phi}{\left( {N + h} \right)}} & {{{- 2}\omega^{e}\cos \quad \phi} - \frac{V^{e}}{N + h}} \\{{{- 2}\omega^{e}\sin \quad \phi} - \frac{2V^{e}\tan \quad \phi}{N + h}} & {- \frac{V^{u}}{M + h}} & \frac{- V_{n}}{M + h} \\{{2\omega^{e}\cos \quad \phi} + \frac{2V^{e}}{N + h}} & \frac{2V^{n}}{M + h} & 0\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad V^{e}} \\{\delta \quad V^{n}} \\{\delta \quad V^{u}}\end{pmatrix}} +}} \\{\quad {R_{b}^{l}\begin{pmatrix}{\delta \quad f_{x}} \\{\delta \quad f_{y}} \\{\delta \quad f_{z}}\end{pmatrix}}}\end{matrix} & (73)\end{matrix}$

[0221] where

[0222] f_(e), f_(n) and f_(u) are three components of the specific forcevector f^(l)=(f_(e) f_(n) f_(u))^(T)=R_(b) ^(l)f^(b)=R_(b) ^(l)(f_(x)f_(y) f_(z))^(T), where f_(x), f_(y) and f_(z) are the specific forcemeasurement outputs of the accelerometers.

[0223] and R is the mean radius of the Earth.

[0224] Measurement Errors

[0225] The sensor measurements can be compensated for the constantbiases in both the FOGs and the accelerometers. However, there are stillsome residual errors which are random in nature. These errors becomesignificant in applications that require long-time surveying. Therefore,these residual random errors are modelled as stochastic processes. Theresidual gyroscopic measurement errors are correlated in time, thereforethey are usually modelled as 1^(st) order Gauss-Markov processes asfollows [Brown et al., 1992; Gelb, 1979]: $\begin{matrix}{\begin{pmatrix}{\delta \quad {\overset{.}{\omega}}_{x}} \\{\delta \quad {\overset{.}{\omega}}_{y}} \\{\delta \quad {\overset{.}{\omega}}_{z}}\end{pmatrix} = {{\left( \quad \begin{matrix}{- \alpha_{x}} & 0 & 0 \\0 & {- \alpha_{y}} & 0 \\0 & 0 & {- \alpha_{z}}\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad \omega_{x}} \\{\delta \quad \omega_{y}} \\{\delta \quad \omega_{z}}\end{pmatrix}} + {\begin{pmatrix}\sqrt{2\alpha_{x}\sigma_{1x}^{2}} \\\sqrt{2\alpha_{y}\sigma_{1y}^{2}} \\\sqrt{2\alpha_{z}\sigma_{1z}^{2}}\end{pmatrix}{W(t)}}}} & (74)\end{matrix}$

[0226] where

[0227] α_(x), α_(y) and α_(z) are the reciprocals of the timecorrelation parameters of the random processes associated with theangular velocities measurements;

[0228] σ_(1x), σ_(1y) and σ_(1z) are the standard deviations of theserandom processes;

[0229] W(t) is unity-variance while Gaussian noise.

[0230] Similarly, the residual accelerometer measurement errors aremodelled as 1^(st) order Gauss-Markov processes as follows:$\begin{matrix}{\begin{pmatrix}{\delta \quad \overset{.}{f_{x}}} \\{\delta \quad \overset{.}{f_{y}}} \\{\delta \quad \overset{.}{f_{z}}}\end{pmatrix} = {{\left( \quad \begin{matrix}{- \beta_{x}} & 0 & 0 \\0 & {- \beta_{y}} & 0 \\0 & 0 & {- \beta_{z}}\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad f_{x}} \\{\delta \quad f_{y}} \\{\delta \quad f_{z}}\end{pmatrix}} + {\begin{pmatrix}\sqrt{2\beta_{x}\sigma_{2x}^{2}} \\\sqrt{2\beta_{y}\sigma_{2y}^{2}} \\\sqrt{2\beta_{z}\sigma_{2z}^{2}}\end{pmatrix}{W(t)}}}} & (75)\end{matrix}$

[0231] where

[0232] β_(z), β_(y) and β_(z) are the reciprocals of the timecorrelation parameters of the random processes associated with thespecific force measurements;

[0233] σ_(2x), σ_(2y) and σ_(2z) are the standard deviations of theserandom processes.

[0234] Equations 71, 72 and 73 for the attitude, position and velocityerrors are combined together with equations 74 and 75 for the FOGs andaccelerometer errors to describe the behavior of the error state e^(l)defined in equation 67. The augmented equation is given as a linearsystem with random system noise as follows:

e ^(.l) =Fe ^(l) +GW(t)  (76)

[0235] where

[0236] F is a 15-by-15 matrix known as the dynamic matrix of the errorstates;

[0237] G is a 15-by-1 vector known as the vector associated with thesystem noise;

[0238] W(t) is unity-variance white Gaussian noise.

[0239] The linear dynamic system shown in equation 75 can be modelled bythe following discrete linear state equation:

e _(k+1) =F _(k+1,k) e _(k) +G _(k) W _(k)  (77)

[0240] where

[0241] e_(k+1) is the error state vector at time t_(k+1);

[0242] e_(k) is the error state vector at time t_(k);

[0243] F_(k+1,k) is the transition matrix of the discrete linear systemdetermined from the dynamic matrix of the error state F asF_(k+1,k)=I+FΔt, where I is the identity matrix and Δt=t_(k+1)−t_(k);

[0244] W_(k) is the random sequence of system noise.

[0245] The problem now is to estimate the error state e_(k) using thediscrete linear state equation given in equation 76. Due to the randomnoise component existing in this linear state equation, its solutionbecomes a filtering problem which is usually addressed using discreteKalman filter theory.

[0246] Aided Strapdown Inertial Navigation Employing Kalman FilteringTheory

[0247] Kalman filtering is a well-known technique for applied optimalestimation especially in the field of inertial navigation and kinematicpositioning [Brown et al., 1992]. Kalman filtering is applied to enhancethe performance of the developed techniques by removing the effect ofresidual random errors during continuous surveying. In addition, Kalmanfilter delivers real-time statistical data related to the accuracy ofthe estimated values. To apply the Kalman filter, an independentmeasurement equation should be specified in addition to the systemequation described in the previous section (equation 76). Thismeasurement equation should be based on an independent measurement ofsuperior accuracy different from the measurements provided by the FOGsand the accelerometers [Hulsing, 1989].

[0248] In MWD borehole surveying, continuous measurement of the lengthof the drill pipe as well as real-time measurement of the BHApenetration rate can respectively provide information about the depth ofthe borehole assembly (i.e. the variable h in the navigation vector) andthe three velocity components V^(e), V^(n) and V^(u) . This independentmeasurement provides what is known as aided strapdown inertialnavigation. FIG. 18 shows a block diagram of the technique of aidedinertial navigation employing Kalman filtering.

[0249] The measurements provided by the FOGs and the accelerometers aresent to the surface of the drilling site with an appropriate telemetrytechnique. These measurements are processed with the surface computersystem to provide the navigation parameters (three position components,three velocity components and three attitude angles). The error statesare then estimated using Kalman filtering which utilizes the measurementof the length of the drill pipe and the BHA penetration rate as anindependent information of superior accuracy. According to the techniqueof aided inertial navigation [Brown and Hwang, 1992], the length of thedrill pipe, L, should be converted first into the corresponding truevertical depth (TVD), h*, as follows:

h*=L cos ζ  (78)

[0250] where

[0251] ζ Is the inclination angle.

[0252] h* Is the true vertical depth as calculated using the length ofthe drill pipe.

[0253] The true vertical depth determined using the length of the drillpipe, h*, is compared to the corresponding value obtained by thecomputational navigation procedure, h. Similarly, the velocitycomponents determined using the measurements of the BHA penetration rateare compares with the corresponding values obtained by the computationalnavigation procedure. The error in between is processed by Kalman filterto provide an estimate of the error states ê which is removed from theoutput of the navigation procedure to deliver corrected inertial output.Therefore, the measurement equation based on the length of the drillpipe as an aiding source is expressed mathematically as follows:

y _(k) =H _(k) e _(k)+ν_(k)  (79)

[0254] where

[0255] Y_(k) is the error in the measurement vector of both the drillpipe length and the penetration rate.

[0256] H_(k) is the design vector of the measurement equation.

[0257] ν_(k) is the discrete measurement random noise.

[0258] Since the design vector, H_(k), relates the measurement error inthe TVD and the velocity components to the error vector, e_(k), it canbe written as $\begin{matrix}{H_{k} = \left( \quad \begin{matrix}0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{matrix}\quad \right)} & (80)\end{matrix}$

[0259] The problem now is to estimate the error states e_(k) based onthe set of independent information y_(k) with the assumption that thesystem noise W_(k) and the measurement noise ν_(k) are uncorrelated. Thecovariance matrix of the system noise Q_(k) and the covariance matrix ofthe measurement noise R_(k) are given as $\begin{matrix}{{E\left\lbrack {W_{k}W_{l}^{T}} \right\rbrack} = \left\{ \begin{matrix}Q_{k} & {k = l} \\0 & {k \neq l}\end{matrix} \right.} & (81) \\{{E\left\lbrack {v_{k}v_{l}^{T}} \right\rbrack} = \left\{ \begin{matrix}R_{k} & {k = l} \\0 & {k \neq l}\end{matrix} \right.} & (82)\end{matrix}$

[0260] with E[W_(k)ν_(l) ^(T)]=0 for all values of k and l.

[0261] The Kalman filter provides a sequential recursive method for theoptimal least mean variance estimation of the system states e_(k). TheKalman filtering is a two step process. The first step is to make aprediction of the error states based on the estimation of the previouserror states. The second step is to update the predicted estimationusing the available measurement from the current step. This approach canbe quantified as [Schwarz et al., 1999]:

[0262] Prediction

ê _(k)(−)=F _(k,k−1) ê _(k−1)(+)

P _(k)(−)=F _(k,k−1) P _(k−1)(+)F _(k,k−1) ^(T) +G _(k−1) Q _(k−1) G_(k−1) ^(T)  (83)

[0263] Where

[0264] ê_(k)(−) is the predicted estimate of the error state at timet_(k);

[0265] ê_(k−1)(+) is the updated estimate of the error state at timet_(k−1);

[0266] P_(k)(−) is the predicted estimate of the covariance matrix ofthe estimation error (P_(k)(−)=E{[ê_(k)(−)−e_(k)][ê_(k)(−)−e_(k)]^(T)}).

[0267] Update

ê _(k)(+)=ê _(k)(−)+K _(k) [y _(k) −H _(k) ê _(k)(−)]

K _(k) =P _(k)(−)H _(k) ^(T) [H _(k) P _(k)(−)H _(k) ^(T) +R _(k)]⁻¹

P _(k)(+)=[I−K _(k) H _(k) ]P _(k)(−)  (84)

[0268] where

[0269] ê_(k)(+) is the updated estimate of the error state at timet_(k);

[0270] K_(k) is the Kalman gain matrix;

[0271] P_(k)(+) is the updated estimate of the covariance matrix of theestimation error (P_(k)(+)=E{[ê_(k)(+)−e_(k)][ê_(k)(+)−e_(k)]^(T)})

[0272] It should be noted that aided inertial navigation based on thelength of the drill pipe as an aiding source will improve the accuracyof the navigation parameters on the vertical channel in a more directway. The parameters of the vertical channel include the true verticaldepth h, the vertical velocity V^(u). On the other hand, the accuracy ofall other navigation parameters (horizontal channel parameters) will notbe significantly affected by the length of the drill pipe as an aidingsource. Alternatively, the independent measurement of the BHApenetration rate will improve the accuracy of the velocity componentsV^(e) and V^(n) and consequently the inclination and roll angles. Inaddition, the navigation parameters of the horizontal channels (East andNorth channels) have another advantage. The errors in the navigationparameters of the horizontal channels are bound due to the Schullereffect and oscillate with a period of 84 minutes (a frequency of about1/5000 known as the Schuller frequency) [Mohammed, 1999; Salychev,1998].

[0273] As mentioned earlier, the technique of aided inertial navigationbased on independent measurement of superior accuracy suits thecontinuous surveying procedure performed inside the radical section ofthe well. During the station-based surveying, the attitude angles aredetermined using the procedure shown on FIG. 17. The values determinedusing this procedure are considered coarse computation and precisecomputation is still needed to improve the estimation accuracy. Thetechnique of precise computation is known as zero velocity update(ZUPT). The ZUPT technique utilizes the condition of zero velocityexisting at each survey station (while the complete surveying setup isstationary) to update the estimate of errors in the Kalman filteringmethod.

[0274] Precise Computation of the Attitude Angles During Station-basedSurveying Employing Kalman Filtering Theory

[0275] The augmented linear state equation (equation 76) describing thebehavior of the error states in the l-frame and its correspondingdiscrete representation is used for station-based precise computation.The only change introduced to suit the station-based surveying is to setthe velocity components V^(e), V^(n) and V^(u) equal to zero. Therefore,the position error state equation is given as follows (see equation 72):$\begin{matrix}{{\delta \quad {\overset{.}{r}}^{l}} = {\left( \quad \begin{matrix}0 & \frac{1}{M + h} & 0 \\\frac{1}{\left( {N + h} \right)\cos \quad \phi} & 0 & 0 \\0 & 0 & 1\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad V^{e}} \\{\delta \quad V^{n}} \\{\delta \quad V^{u}}\end{pmatrix}}} & (85)\end{matrix}$

[0276] Similarly, the velocity error state equation and the attitudeerror state equation are given respectively as follows (see equation 73and 71): $\begin{matrix}\begin{matrix}{\begin{pmatrix}{\delta \quad {\overset{.}{V}}^{e}} \\{\delta \quad {\overset{.}{V}}^{n}} \\{\delta \quad {\overset{.}{V}}^{u}}\end{pmatrix} = \quad {{\left( \quad \begin{matrix}0 & 0 & 0 \\0 & 0 & 0 \\0 & 0 & \frac{2\quad \gamma}{R}\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad \phi} \\{\delta \quad \lambda} \\{\delta \quad h}\end{pmatrix}} +}} \\{\quad {{\left( \quad \begin{matrix}0 & {2\omega^{e}\sin \quad \phi} & {{- 2}\omega^{e}\cos \quad \phi} \\{{- 2}\omega^{e}\sin \quad \phi} & 0 & 0 \\{2\omega^{e}\cos \quad \phi} & 0 & 0\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad V^{e}} \\{\delta \quad V^{n}} \\{\delta \quad V^{u}}\end{pmatrix}} +}} \\{\quad {{\left( \quad \begin{matrix}0 & f_{u} & {- f_{n}} \\{- f_{u}} & 0 & f_{e} \\f_{n} & {- f_{e}} & 0\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad \theta} \\{\delta \quad \varphi} \\{\delta \quad \psi}\end{pmatrix}} + {R_{b}^{l}\begin{pmatrix}{\delta \quad f_{x}} \\{\delta \quad f_{y}} \\{\delta \quad f_{z}}\end{pmatrix}}}}\end{matrix} & (86) \\\begin{matrix}{\begin{pmatrix}{\delta \quad \overset{.}{\theta}} \\{\delta \quad \overset{.}{\varphi}} \\{\delta \quad \overset{.}{\psi}}\end{pmatrix} = \quad {{\left( \quad \begin{matrix}0 & 0 & 0 \\{\omega^{e}\sin \quad \phi} & 0 & 0 \\{{- \omega^{e}}\cos \quad \phi} & 0 & 0\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad \phi} \\{\delta \quad \lambda} \\{\delta \quad h}\end{pmatrix}} +}} \\{\quad {{\left( \quad \begin{matrix}0 & \frac{1}{M + h} & 0 \\\frac{- 1}{N + h} & 0 & 0 \\{{- \tan}\quad \phi} & 0 & 0\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad V^{e}} \\{\delta \quad V^{n}} \\{\delta \quad V^{u}}\end{pmatrix}} +}} \\{\quad {{\left( \quad \begin{matrix}0 & {\omega^{e}\sin \quad \phi} & {{- \omega^{e}}\cos \quad \phi} \\{{- \omega^{e}}\sin \quad \phi} & 0 & 0 \\{\omega^{e}\cos \quad \phi} & 0 & 0\end{matrix}\quad \right)\begin{pmatrix}{\delta \quad \theta} \\{\delta \quad \varphi} \\{\delta \quad \psi}\end{pmatrix}} + {R_{b}^{l}d}}}\end{matrix} & (87)\end{matrix}$

[0277] The measurement error states of both the accelerometer and theFOGs are described as first order Gauss-Markov process and they areexpressed as follows: $\begin{matrix}{\begin{pmatrix}{\delta {\overset{.}{\omega}}_{x}} \\{\delta {\overset{.}{\omega}}_{y}} \\{\delta {\overset{.}{\omega}}_{z}}\end{pmatrix} = {{\begin{pmatrix}{- \alpha_{x}} & 0 & 0 \\0 & {- \alpha_{y}} & 0 \\0 & 0 & {- \alpha_{z}}\end{pmatrix}\begin{pmatrix}{\delta\omega}_{x} \\{\delta\omega}_{y} \\{\delta\omega}_{z}\end{pmatrix}} + {\begin{pmatrix}\begin{matrix}\sqrt{2\alpha_{x}\sigma_{1x}^{2}} \\\sqrt{2\alpha_{y}\sigma_{1y}^{2}}\end{matrix} \\\sqrt{2\alpha_{z}\sigma_{1z}^{2}}\end{pmatrix}{W(t)}\quad {and}}}} & (88) \\{\begin{pmatrix}{\delta {\overset{.}{f}}_{x}} \\{\delta {\overset{.}{f}}_{y}} \\{\delta {\overset{.}{f}}_{z}}\end{pmatrix} = {{\begin{pmatrix}{- \beta_{x}} & 0 & 0 \\0 & {- \beta_{y}} & 0 \\0 & 0 & {- \beta_{z}}\end{pmatrix}\begin{pmatrix}{\delta \quad f_{x}} \\{\delta \quad f_{y}} \\{\delta \quad f_{z}}\end{pmatrix}} + {\begin{pmatrix}\begin{matrix}\sqrt{2\beta_{x}\sigma_{2x}^{2}} \\\sqrt{2\beta_{y}\sigma_{2y}^{2}}\end{matrix} \\\sqrt{2\beta_{z}\sigma_{2z}^{2}}\end{pmatrix}{W(t)}}}} & (89)\end{matrix}$

[0278] The augmented error linear state equation during thestation-based surveying is obtained by combining equations (85-89). Thisaugmented equation is similar to equation 77 but with differentdefinitions of the dynamic matrix F and the system noise vector G whichare given as: $\begin{matrix}{F = \left( \quad \begin{matrix}0 & 0 & 0 & 0 & \frac{1}{M + h} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & \frac{1}{\left( {N + h} \right)\cos \quad \phi} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & {2\omega^{e}\sin \quad \phi} & {{- 2}\omega^{e}\cos \quad \phi} & 0 & f_{u} & {- f_{n}} & 0 & 0 & 0 & r_{11} & r_{12} & r_{13} \\0 & 0 & 0 & {{- 2}\omega^{e}\sin \quad \phi} & 0 & 0 & {- f_{u}} & 0 & f_{e} & 0 & 0 & 0 & r_{21} & r_{22} & r_{23} \\0 & 0 & \frac{2\gamma}{R} & {2\omega^{e}\cos \quad \phi} & 0 & 0 & f_{n} & {- f_{e}} & 0 & 0 & 0 & 0 & r_{31} & r_{32} & r_{33} \\0 & 0 & 0 & 0 & \frac{1}{M + h} & 0 & 0 & {\omega^{e}\sin \quad \phi} & {{- \omega^{e}}\cos \quad \phi} & r_{11} & r_{12} & r_{13} & 0 & 0 & 0 \\{\omega^{e}\sin \quad \phi} & 0 & 0 & \frac{- 1}{\left( {N + h} \right)} & 0 & 0 & {{- \omega^{e}}\sin \quad \phi} & 0 & 0 & r_{21} & r_{22} & r_{23} & 0 & 0 & 0 \\{{- \omega^{e}}\cos \quad \phi} & 0 & 0 & \frac{{- \tan}\quad \phi}{\left( {N + h} \right)} & 0 & 0 & {\omega^{e}\cos \quad \phi} & 0 & 0 & r_{31} & r_{32} & r_{33} & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- \alpha_{x}} & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- \alpha_{y}} & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- \alpha_{z}} & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- \beta_{x}} & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- \beta_{y}} & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- \beta_{z}}\end{matrix}\quad \right)} & (90)\end{matrix}$

[0279] and

G=(0 0 0 0 0 0 0 0 0 {square root}{square root over (2α_(x)σ_(1x) ²)}{square root}{square root over (2α_(y)σ_(1y) ²)} {square root}{squareroot over (2α_(z)σ_(1z) ²)} {square root}{square root over (2β_(x)σ_(2x)²)} {square root}{square root over (2β_(y)σ_(2y) ²)} {squareroot}{square root over (2β_(z)σ_(2z) ²)})^(T)  (91)

[0280] In the case of station-based surveying, the measurement equation,y_(k)=H_(k)e_(k)+ν_(k), is completely different from the aided inertialnavigation technique discussed in the previous section. There is noaiding source that can provide a measurement of superior accuracy.Alternatively, we employ the condition of zero velocity to construct themeasurement equation. This condition is introduced because the bottomhole assembly is at a complete rest at each survey station. In thesecircumstances the error measurement vector y_(k) is given as follows

y _(k) =V ^(l) −V ⁰  (92)

[0281] where

[0282] V^(l) is the velocity vector expressed in the l-frame,V^(l)=(V^(e) V^(n) V^(u))^(T);

[0283] V⁰ is the zero velocity vector, V⁰=(0 0 0) ^(T).

[0284] Since the error measurement vector, y_(k), expresses an error invelocity, the measurement equation is given similarly to equation 79 butwith the following definition of the design matrix H, $\begin{matrix}{H = \left( \quad \begin{matrix}0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{matrix}\quad \right)} & (93)\end{matrix}$

[0285] Kalman filtering described in equations 83 and 84 can then beapplied but using the new definitions for F and H (the dynamic and thedesign matrices respectively), as well as the new definitions of the Gand y_(k) (the system noise and error measurement vectors).

[0286] The precise computation technique described in this sectionimproves the estimation accuracy by removing the residual random noisecomponents existing in the accelerometer and the FOG measurements. Sincethe ZUPT technique is directly related to the velocity components of thebottom hole assembly, the velocity errors, δV^(e), δV^(n) and δV^(u),will be significantly affected in a more direct way than the positionerrors, δφ, δλ and δh, and the attitude errors, δθ, δφ and δΨ. Theprecise computation usually takes about 10 to 15 minutes for the Kalmanfilter parameters to converge. The exact time needed for convergencedepends mainly on the sensor accuracy and the level of the random noiseat their outputs.

[0287] As will be apparent to those skilled in the art, variousmodifications, adaptations and variations of the foregoing specificdisclosure can be made without departing from the scope of the inventionclaimed herein.

REFERENCES

[0288] The contents of the following references are hereby incorporatedherein as if reproduced herein in their entirety.

[0289] Britting K. R.: “Inertial navigation system analysis;” John Wileyand Sons, New York, 1971.

[0290] Brown R. G. and Hwang P. Y. C.: “Introduction to random signals;”John Wiley and Sons, New York, 1992.

[0291] Conti P. F.: “Controlled horizontal drilling;” Proceedings ofSPE/IADC 1989 Drilling Conference, New Orleans, La. pp:749-754, Feb.28-Mar. 3, 1989.

[0292] Disperio R.: “Method of correcting for axial error components inmagnetometer reading during wellbore survey operation; ” U.S. Pat. No.5,452,518,1995.

[0293] Engebretson H.: “Error reduction in compensation of drill stringinterference for magnetic survey tools;” U.S. Pat. No. 5,155,916,1992.

[0294] Gelb A.: “Applied optimal estimation;” MIT Press, Cambridge,England, 1989.

[0295] Grindord S. and Wolf C.: “Calculation of NMDC length required forvarious latitudes developed from field measurements of drill stringmagnetization;” Proceedings of SPE/IADC 1983 Drilling Conference, NewOrleans, Tex. pp: 217-224, Feb. 20-23, 1983.

[0296] Hartman R.: “Method of qualifying a borehole survey:” U.S. Pat.No. 5,787,997, 1998.

[0297] Helm W.: “Method and apparatus for measurement of azimuth of aborehole while drilling;” U.S. Pat. No. 5, 012,412, 1991.

[0298] Hulsing R.: “Borehole survey system utilizing strapdown inertialnavigation;” U.S. Pat. No. 4,812,977, 1989.

[0299] Lefevre H.: “The fiber optic gyroscope;” Artech-House, (USA),1993.

[0300] Kim B. Y.: “Optical fiber rotation sensing;” Academic press Inc,New York (USA) 1994.

[0301] Merhav, S.: “Aerospace sensor systems and applications;”Springer, N.Y. (USA), 1993.

[0302] Mohamed A. H.: “Optimizing the estimation procedure in INS/GPSintegration for kinematic applications;” Ph.D. thesis, Department ofGeomatics Engineering, University of Calgary, UCGE report # 20127, April1999.

[0303] Molnar, D.O. “Borehole Inertial Guidance System” U.S. Pat. No.4,542,647.

[0304] Nichelson J.: “Method for determining borehole direction;” U.S.Pat. No. 5,435,969, 1995.

[0305] Noureldin A., Mintchev M., Irvine-Halliday D. and Tabler H.:“Computer modelling of microelectronic closed loop fiber opticgyroscope;” Proceedings of the IEEE Canadian Conference on Electricaland Computer Engineering; Edmonton, Canada, pp: 633-638, May 9-12, 1999.

[0306] Noureldin A., Tabler H., Irvine-Halliday D. and Nintchev M.:“Quantitative study of the applicability of fiber optic gyroscopes formeasurement-while-drilling borehole surveying;” Journal of the Societyof Petroleum Engineers (SPE), accepted for publication, May 2000.

[0307] Parkinson W.D.: “Introduction to Geomagnetism;” Scottish AcademicPress, 1983.

[0308] Rehm W. A. Garcia A. and Cia S. A.: “Horizontal drilling inmature oil fields;” Preceedings of SPE/IADC 1989 Drilling Conference,New Orleans, La., pp: 755-764, Feb. 29-Mar. 3, 1989.

[0309] Russel, M. K. and Russel A. W.: “Surveying of boreholes;” U.S.Pat. No. 4,163,324, 1979.

[0310] Salychev O.: “Intertial systems in navigation and geophysics;”Bauman MSTU Press, Moscow, 1998.

[0311] Schwarz K. P. and Wei M.: “INS/GPS integration for geodeticapplications;” Lecture notes for ENGO 623, Department of GeomaticsEngineering at the University of Calgary, 1999.

[0312] Schwarz K. P. and Wei M.: “A framework for modelling kinematicmeasurements in gravity field applications;” Bulleting Geodesique, v 64(4), 331-346, October 1990.

[0313] Sheills G., Kerridge D.: “Borehole surveying;” U.S. Pat. No.6,021,577, 2000.

[0314] Skillingstad T.: “At-bit inclination measurements improvesdirectional drilling efficiency and control;” Proceedings of SPE/IADC2000 Drilling Conference, New Orleans, La., Feb. 23-25, 2000.

[0315] Stephenson M. A. and Wilson H.: “Improving quality control ofdirectional survey data with continuous inertial navigation;” SPEDrilling Engineering, v 7 (2), pp: 100-105, June 1992.

[0316] Thorogood J. L.: “How to specify and implement well surveys;”World Oil, July 1986.

[0317] Thorogood J. L. and Knott, D. R.: “Surveying techniques with asolid state magnetic multi-shot device;” Proceedings of SPE/IADC 1989Drilling Conference, New Orleans, La., pp: 841-856, Feb. 28-Mar. 3,1989.

[0318] Thorogood J. L.: “Instrument performance models and theirapplication to directional surveying operations;” SPE DrillingEngineering, v 5 (4), pp: 294-298, December 1990.

[0319] Titterton D. H. and Weston, J. L.: “Strapdown inertial navigationtechnology;” Peter Peregrinus Ltd., London, UK, 1997.

[0320] Trowsdale L.: “Borehole survey method and apparatus for drillingsubstantially horizontal boreholes;” U.S. Pat. No. 4,361,192, 1982.

[0321] Van Dongen J., Maekiaho L.: “Method for determining the azimuthof a borehole;” U.S. Pat. No. 4,682,421, 1987.

[0322] Yakowitz S. and Szidarovszky F.: “An introduction to numericalcomputations;” 2^(nd) edition, Macmillan publication company, New York,1989.

What is claimed is:
 1. A continuous measurement-while-drilling surveyingapparatus for surveying the drilling progress of a bottom hole assembly(“BHA”) having a tool-spin axis, said apparatus comprising: (a) a firstfiber-optic gyroscope for generating a first angular rotation signalrepresenting angular rotation of the BHA about the tool-spin axis; (b) asecond fiber-optic gyroscope mounted within the housing for generating asecond angular rotation signal representing angular rotation of the BHAabout a second axis normal to the tool-spin axis; (c) accelerometermeans for generating three acceleration signals representing thecomponents of acceleration of the BHA along three mutually orthogonalaxes; (d) first processing means responsive to the acceleration signalsfor determining the angle of the BHA away from the vertical and forgenerating a third angular rotation signal representing rotation of theBHA about an axis normal to the sensitive axes of the first and secondgyroscopes; (e) second processing means responsive to the first, secondand third angular rotation signals and the acceleration signals fortransforming signals representing movement of the BHA in a BHAcoordinate system to a earth local-level coordinate system; and (f)third processing means operatively connected to the second processingmeans for determining the orientation of the BHA, determining thevelocity changes of the BHA, updating the velocity components of the BHAand updating the position components of the BHA.
 2. The apparatus ofclaim 1 wherein the first gyroscope is toroidal.
 3. The apparatus ofclaim 1 wherein the second axis is aligned with a forward direction ofthe BHA.
 4. The apparatus of claim 1 further comprising a Kalman filteroperatively connected to the third processing means.
 5. A continuousmeasurement-while-drilling surveying apparatus for surveying thedrilling progress of a bottom hole assembly (“BHA”) having a tool-spinaxis, said apparatus comprising: (a) a toroidal fiber-optic gyroscopesecured within or adjacent to the BHA such that the gyroscope defines acentral passageway for the flow of drilling fluid; and (b) processingmeans for receiving the output of the gyroscope and producing a signalrepresentative of the angular velocity of the BHA about the tool-spinaxis.
 6. A method of continuous MWD surveying of a wellbore whichincludes, on completion, a vertical section, a radical section and ahorizontal section, using a BHA comprising a first gyroscope having itssensitive axis aligned with the tool spin axis, a second gyroscopehaving its sensitive axis on an axis normal to the tool spin axis andthree accelerometers arranged in three mutually orthogonal directions,said method comprising the steps of: (a) receiving the rotation signalsfrom the first gyroscope and the second gyroscope and the accelerationsignals from the three accelerometers; (b) establishing the desiredazimuth direction; (c) determining the pitch angle from theaccelerometer signals; (d) determining the rotation rate about an axisnormal to the sensitive axes of the two gyroscopes by determining thetime rate of change of the pitch angle; (e) determining the effect ofthe earth rotation and local-level (l-frame) change of orientation tothe angular changes monitored along the three axes and removing thesetwo effects from the rotation rate signals; (f) transforming therotation signals from a BHA coordinate frame to the earth local-levelcoordinate frame; (g) determining the pitch, roll and azimuth of theBHA; (h) determining the time rate of change of the velocity of the BHAexpressed in the l-frame and update the velocity components of the BHA;and (i) updating the position components of the BHA.
 7. The method ofclaim 6 wherein the rotation signal from the first gyroscope is used tomonitor changes in the azimuth of the BHA during drilling of thevertical section and an initial portion of the radical section of thewell and is thereafter used to monitor changes in the roll angle of theBHA.
 8. The method of claim 6 wherein the rotation signal from thesecond gyroscope is used to monitor changes in the roll angle of the BHAduring drilling of the vertical section and an initial portion of theradical section of the well and is thereafter used to monitor changes inthe azimuth of the BHA.
 9. The method of claim 6 further comprising thestep of removing bias terms from the acceleration signals beforedetermining the pitch angle.
 10. The method of claim 6 furthercomprising the step of removing gyroscope drift from the rotationsignals.
 11. The method of claim 6 further comprising the step ofcorrecting inertial output of the method by estimating error statesusing Kalman filtering and independent information as an aiding source.12. The method of claim 11 wherein the independent information comprisesmeasurements of the length of drill pipe and BHA penetration rate. 13.The method of claim 7 wherein the initial portion of the radical sectionof the well is that portion having an inclination less than about 45°.14. The method of claim 13 wherein the inclination is less than about30°.
 15. The method of claim 14 wherein the inclination is less thanabout 20°.
 16. A method of continuous MWD surveying the progress of abottom hole assembly having a tool-spin axis, an azimuth angle and atool face angle, said method comprising the steps of: (a) providing atoroidal fiber-optic gyroscope within the BHA, wherein said gyroscopeisoriented horizontally when the BHA is oriented vertically and defines acentral passage for the flow of drilling fluid; and (b) monitoring therotation of the BHA about the tool-spin axis with the gyroscope.
 17. Themethod of claim 15 wherein the rotation of the BHA about the tool-spinaxis is used to determine the azimuth angle of the BHA during a verticalsection of the wellbore and an initial portion of a radical section ofthe wellbore.
 18. The method of claim 16 wherein the rotation of the BHAabout the tool-spin axis is used to determine the tool face anglesubsequent to the initial portion of the radical section of thewellbore.